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Tables are presented of thermodynamic properties of ethylene in the range 279-300 K in 
temperature, 5.75 - 10.5 mol/dm? in density, which range includes the critical point. The 
tables presented here are based on the critical-point scaling laws and incorporate the 
critical anomalies in a proper fashion. The tables thus complement the formulation of the 
equation of state of fluid ethylene by McCarty and Jacobsen (NBS Tech. Note 1045, 1981) 
which does not claim accuracy near the critical point. The predictions of the present 
formulation are compared with four sets of recent PVT data, and with speed-of-sound data. 
Tables are presented of pressure, energy, enthalpy, entropy, specific heats and speed of 
sound as function of temperature along finely-spaced isochores. The computer program 
required for table generation is included. Even if the surface were perfect, the reliabil- 
ity of densities calculated at experimental pressures and temperatures of limited accuracy 
declines rapidly as the critical point is approached. Contour plots in P-T space are 
presented of regions to be avoided in custody transfer for given uncertainties in pressure, 
temperature and sample composition. 


Key words: critical region; custody transfer; density; enthalpy; ethylene; 
equation of state; impurity; scaling laws; specific heat; 
speed of sound; tabular values; thermodynamic properties. 
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te Introduction 


A formulation is presented of the thermodynamic behavior of ethylene in a region around 
the critical point extending from 5.75 to 10.5 mol/dm? in density and from 279 to 300 K in 
temperature. A correct formulation of this near-critical region is more urgent for ethylene 
than for most other fluids, because in the climatic conditions of the United States the 
fluid transported in pipes can assume states located in this region. It is well recognized 
that most engineering equations of state are deficient in representing near-critical states 
[1,2]. Our formulation of the thermodynamic behavior of ethylene derives from a thermodyna- 
mic potential which obeys the critical-point scaling laws; it contains a first correction to 
scaling and also a revision to scaling designed to accommodate the lack of symmetry of the 
gas-liquid phase transition. The adjustable parameters in the potential were obtained by a 
fit to the PVT data of Hastings et al. [3] and to the speed-of-sound data of Gammon [4]. 

The advantage of our formulation is that the critical anomalies that are common to all 
fluids are properly incorporated. For instance, the flatness of the critical isotherm, 
which is hard to describe by analytic equations of state, is built into our potential. Its 
strength is particularly striking when it comes to representing the specific heat at con- 
stant volume, Cys and the speed of sound, C.. The weak divergence of C, cannot be described 
by analytic free energy surfaces. This divergence, however, causes the speed of sound to 
approach the value zero at the critical point, a feature properly incorporated in our poten- 
tial. The excellent speed-of-sound data of Gammon therefore form a crucial test for the 
adequacy of our thermodynamic surface. We stress that, by the use of so-called parametric 
variables, we have closed-form expressions for all thermodynamic properties. 

We will describe how our potential is constructed, and proceed to compare it with 
recent experimental PVT data from three sources and with vapor pressures and coexisting 
densities from two sources. Comparisons are also made with the predictions of the global 
equation of state of ethylene recently published by McCarty and Jacobsen [5]. From these 
comparisons it will be evident that our potential will serve as a useful and necessary 
supplement in a region where the global equation is essentially deficient. In appendices, 
we present the equations defining our potential and a comparison of all data points me 
tioned aboye with the predictions of our surface. A table of the thermodynamic properties 
pressure, energy, entropy, enthalpy, specific heats Cp and Cy and speed of sound Co; is 
presented along finely-spaced isochores with temperature as an entry. A table of saturation 
properties is included. The computer programs for generation of these tables are included. 
They contain an inversion routine which permits transformation from pressure to density as 
an independent variable. 

For the convenience of the user concerned with custody transfer, we conclude the main 
body of the paper with a discussion of the uncertainties introduced in the density calcu- 
lated from our surface due to limited precision in the measurement of pressure and tempera- 
ture, and due to the presence of impurity. Contour plots in P-T space indicating regions in 
PT space to be avoided if a certain precision in density is desired, are included in this 
paper. Our work has been part of a joint industry-U.S. Government project to determine the 
thermodynamic properties of ethylene. 


2. Fundamental Equation 


At a critical point, thermodynamic behavior is anomalous. The critical point is char- 
acterized by the facts that the compressibility Kaos es (1/V) (o¥/ oP) + and the specific heat 


Cy diverge and the coexisting densities Po Py becomes equal. We write, asymptotically near 
the critical point 


eS ae es 


Ao* = BlaT*|® (1) 


Here asterisks denote reduced properties, while AT* = CT i/Te and Ap* = (o,-0,)/20,. The 
critical exponent y equals 1, g equals 1/2 for classical equations of state of the type of 
Van der Waals' equation. In real fluids, y equals about 1.24 and g is of the order 1/3. 
The critical-point scaling laws are constructed so that they yield the correct critical 
exponents; although one can define as many critical exponents as one sees fit, only two of 
them are independent. Their values are given by theory and they are the same for all fluids. 

Our fundamental equation incorporates the scaling laws. The thermodynamic potential 
used in our work is that introduced by Ley-Koo and Green [6]. It is the dependent variable 
p = P*/T* as a function of » = p*/T* and T = -]/T*. Here P is the pressure, »p is the chemi- 
cal potential and T the absolute temperature. Asterisks denote reduced properties: P* = P/P.> 
T* = We and u* = up /P os where Pe is the critical pressure, ile the critical temperature 
and Po the critical density. In terms of this potential, we have 

(aP/9T) ae Wie, 


(aP/au)z= 6 (2) 


where U is the reduced energy density U = U*/V* and p the reduced mass density, while the 
other thermodynamic functions: Helmholtz energy density A = A*/V*, enthalpy density 
H = H*/V* and entropy density S = S*/V* follow from 
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Since all variables are intensive, the coexistent phases collapse onto a curve, which, in 
the fi, T plane, can be parametrized by the analytic relation 


2 aig? = Wig 
HC) Sts eB uy AT, 


where AT = (eT MMe In scaling, the physical variables u, T are transformed to a new 


coordinate system, one axis of which is the curve uy (1); while the other intersects the 
first at the critical point. The two scaling variables are linear combinations of nau and 
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where a is a constant related to the scale of up. The term cAy breaks the liquid-vapor 
symmetry. 

The potential P is decomposed into two parts, one part, ae that is regular in the 
scaling variables and one, AP, that contains the critical anomalies with correct power-law 


dependence. We have assumed for Pe : 


g 
Bish = Pa(T) fhe Pan AM IE 
~ ~ 4 (5) 
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AP, in scaled form, is given by 
u u 
Fi 2-o yu 2-atAy u 
AP = ak |u,| ad aaah ak, |u, | ( 2 (6) 
lu, | ju, 


The first term represents asymptotic scaling, the second a correction to scaling. The 
critical exponents a,8,6 are those characterizing the divergence of the specific heat at 
constant volume Cy» the shape of the coexistence curve, and the shape of the critical isotherm, 
respectively. The leading term contains one more adjustable constant koe while the function 
Io is universal, i.e., the same for all fluids. The correction term contains a new inde- 
pendent critical exponent Ay» also given by theory, and a new free constant k,. The 
function g, is the same for all fluids. 

All thermodynamic functions can be given in algebraically closed form by the device of 
transformation to parametric variables [7] 


Ap = au = ao (eae). 


u, =r (1-b202) (7) 
The variable r measures "distance from the critical point", while 6 is measured along a 
contour of constant r. 6 reaches the values +1 on the phase boundary; b* is a constant 
predicted to be universal. The critical anomalies are all imbedded in the r-dependence 
while the dependence of the thermodynamic functions on 6 is always in the form of simple 
polynomials or their ratios. All relevant equations are given in Appendix A. Note that the 
functions are listed separately for the one- and two-phase regions. 

The constants defining our potential were obtained as follows. The critical exponent 
values a,8,7,6,A, are those predicted by theory [8]. The critical parameters, the five 
parameters in the scaling function: a, Ko k,, C, b* and the four "background" parameters 


Pre Das De and Ban were determined by a fit to the PVT data of Hastings et al. [3], while we 
imposed the condition that the slope of the vapor pressure curve at the critical point be 
that of the vapor pressure data of Levelt Sengers and Hastings [9], that the critical point 
lie on the vapor pressure curve and that the diameter of the coexistence curve be that 
measured by Douslin and Harrison [10]. This implies that the values of ie and Re are con- 
strained to be on the vapor pressure curve and cannot be fitted independently. Two back- 
ground parameters jiz,ji3 were determined by a fit to the speed-of-sound data of Gammon [4]. 
The other two, Ue and j1,, are related to the zeros for energy, entropy and enthalpy, and 
were fixed by tying the surface to that of McCarty and Jacobsen [5] at one point, 298.15 K 
and 6 mol/dm?. The parameters characterizing our potential are listed in Table I. The 
range of validity of the potential is indicated in Fig. 1. 


3. The Equation of State 


We compare the predictions of our potential for the equation of state of ethylene with 
experiment. The experimental data are those of Hastings et al. [3] for a 99.993% pure 
sample, to which the potential was fitted, those of Douslin and Harrison [10], Trappeniers 
et al. [11] and Thomas and Zander [12,13]. Older data have been adequately reviewed in the 
IUPAC formulation of Angus et al. [14], and in the reviews of Date [15] and Vashchenko [16]. 

The potential represents the data of Hastings et al. with a reduced chi-square of 1.3 
if the uncertainties are estimated as 5 x 107° MPa in pressure, 3 x 10-2 mol/dm? in density 
and 5 x 107* K in temperature. It is to be kept in mind that our data were fitted to these 
close tolerances even though the deviation plots of Figs. 2-13 show occasionally large 
density deyiations. The percentage differences in density between the data and our surface, 
assuming the pressures to be the same, and, alternatively, the percentage differences in 
pressure, assuming the densities to be the same, are shown for the various isotherms in 
Figs. 2-13. At 298.15K (Fig. 2) the surface represents our own data to well within 0.1% in 
density. The agreement between the data sets of [3], [10] and [11] is almost within the 
0.1% tolerance originally set as a goal for the U.S. Government - Industry joint project on 
properties of ethylene. The data set of Thomas and Zander [12,13] shows considerably more 
Scatter than the other three sets, without, however, showing systematic departures from our 
formulation. The pressure deviations are similarly small at this temperature, with a ten- 
dency to grow in the highly compressed fluid, as is to be anticipated [Fig. 3]. When the 
critical point is approached, the density departures between the data sets, and their 
departures from the surface, increase dramatically, while the pressure deviations remain 
more or less the same (Figs. 4-10). At the temperatures from 0.7-3.7 K above the critical 
temperature, the density departures reach values larger than 0.5% and on isotherms spanning 
the critical the largest departure is 15%. The pressure deviations, however, remain as 
well-behaved as at 298.15K if not better, the departures being of the order of 0.02% in the 
range around critical (Figs. 11-13). This is a good illustration of the difference between 
"density" and "field" variables in the vicinity of a critical point. ("Density" variables 
are those that are different in coexisting phases while "field" variables are the same.) 
Due to the strong divergences of the compressibility and expansion coefficients, smal] 


systematic differences in temperature and pressure scales between the various laboratories 
are amplified when the densities are compared. Furthermore, the effect of impurity is 
largest at the critical point, and small differences in sample composition can lead to 
large variations in density. 

A study of the plots Figs. 2-13 will convince the reader that most of the region 
covered in this report should be avoided if an accuracy of 0.1% in density is desired in 
custody transfer. Even if our surface were perfect, the finite precision of industrial 
pressure and temperature measurement will result in large errors if the density is calcu- 
lated from the equation of state. 

Appendix B, Table II, contains the deviations of the PVT data sets [3], [10], [11], 
[12], [13] and [14] from the present formulation. 


4, Coexistence Properties 


Vapor pressures from 279.15 K to the critical temperature are compared in Fig. 14. 
Data sources are those of Levelt Sengers and Hastings [9] and those of Douslin and Harrison 
[10]. These two data sets were shown to agree within combined error of 107% MPa at tempera- 
tures below 240 K [9]; the departures are systematic at higher temperatures, and reach a 
maximum of 7 x 107% MPa at the critical point. A satisfactory explanation for the differ- 
ence has not been given. It is larger than expected on the basis of the agreement of pres- 
sure scales (1 in 20,000) and temperature scales (1 mK) and not quite of the form expected 
for a difference in purity [7]. Nevertheless, the agreement between the data sets [9] and 
[10], better than 2 parts in 10', is an order better than that between these sets and the 
older data (cf. [10], p. 311). Table III in Appendix B contains the deviations of the vapor 
pressures data sets [9] and [10] from our formulation. 

Coexisting densities of Refs. [3] and [10] are compared with our formulation in 
Fig. 15. The agreement between the two data sets and the formulation, to +0.02 mol/dm?, is 
satisfactory in this region close to the top of the coexistence dome. 


De Speed of Sound 


As we mentioned before, we determined the parameters Ho» ls in our potential from a fit 
to Gammon's speed-of-sound data [4], accepting his density and temperature values as inde- 
pendent variables. It should be stressed that the anomalous contribution to the speed of 
sound that drives speed of sound to zero at the critical point has. been completely predicted 
from the part of the potential that was determined by the PVT data alone. A substantial loss 
of accuracy is to be expected when second derivatives are calculated from the potential. We 
therefore consider the result of our procedure, a prediction of the speed of sound that is 
generally well within 1% from experiment, as very satisfactory. A global view of the speed 
of sound within the range of our correlation is shown in Fig. 16. Fig. 17 shows the depar- 
tures of the speed of sound from our surface. It should be mentioned that we had to omit a 
few of Gammon's data points because our surface predicted them to be inside the two-phase 
region. Also, we omitted the data points on the critical isotherm within + 10% from the 
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critical density, these points being heavily affected by dispersion. 

Fig. 16 and 17 reveal one defect of our model: on subcritical isotherms, the speed of 
sound, instead of continuing to fall as the phase boundary is approached, passes through a 
minimum and then experiences a minute increase in the last percent of density from coexis- 
tence. We believe this behavior to be unphysical. It was not present when we fitted our 
potential to the experimental data for steam. We have, however, not been able to eliminate 
this artifact by the device of clamping more tightly on the PVT data near the phase boundary, 
nor by reducing the range of the fit. It is likely that speed-of-sound data of the quality 
of [4] may contain more information than even high-quality PVT data; it would be interesting 
to use the sound data to determine some of the parameters in the scaled part of the poten- 
tial. 

Table IV in Appendix B contains a listing of the departures of the experimental data 
from the present formulation. 


6. Other Formulations 


A recent comprehensive correlation of the properties of ethylene, in the form of an 
equation of state based on the same data sets considered in this paper, is that of McCarty 
and Jacobsen [5]. It was preceded by several other correlations, two of which we will 
consider here: the IUPAC formulation by Angus et al. [14], which formulation is based on a 
part-graphical, part-analytic representation of the PVT properties of ethylene obtained 
prior to 1972; and the correlation obtained by graphical means by Harrison and Douslin [7] 
of their own data. The solid curves in Figs. 2-13 indicate the departures of the McCarty - 
Jacobsen formulation from our own. Since the Helmholtz function of McCarty and Jacobsen is 
analytic, it has difficulty representing the flatness of near-critical isotherms, as was 
recognized by the authors. At 298.15 K the departures, though systematic, are limited to a 
maximum of 0.3% in density, but as the temperature approaches the critical the departures 
become more severe. At temperatures within 1 K from critical where our formulation agrees 
with all data to well within 0.1% in pressure, the McCarty-Jacobsen formulation departs by 
several tenths of a percent in pressure in a systematic fashion. These differences are 
magnified when derived properties are considered. An example is shown in Fig. 17 for the 
enthalpy of ethylene, where enthalpies according to three formulations are compared along 
the 283.15 K, 288.15 K, and 298.15 K isotherms. The graphical method of Harrison and 
Douslin, if applied with care to excellent data, should be expected to give better results 
than can be obtained with a fit to an empirical thermodynamic surface that is seldom quite 
true to the data. A possible exception is the near-critical regime: The anomalous curva- 
ture of isochores may be missed by graphical techniques and by analytic surface alike, but 
is correctly incorporated in our potential. The intercomparison in Fig. 18 shows that 
the Harrison and Douslin enthalpies show no density-dependent departures from our surface 
and haye, therefore, been very carefully evaluated. At 298.15 K the offset is about 8 J/mol 
and could have been eliminated had we tied our enthalpy at 298.15 K, 6 mol/dm? to that of 
Harrison and Douslin, rather than to that of McCarty and Jacobsen. At the other two tem- 
peratures a systematic offset of 20 J/mol is observed. We do not know why a 12 J/mol 
enthalpy difference is picked up between 298.15 and 288.15 K. 
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At 298.15 K our enthalpies agree closely with those of McCarty and Jacobsen, but the 
density dependence of the two formulations is different. Increasing departures build up as 
the critical point is approached. 

The departures of the IUPAC formulations show strong systematic trends. This formulation 
was based on older data and was never claimed to possess high quality in the critical region. 

First and second derivatives according to the McCarty - Jacobsen formulation are com- 
pared with our own in Figs. 19, 20. The first derivatives U, H and S have somewhat dif- 
ferent density dependence in the two formulations. Moreover, appreciable departures build up 
around the critical density and near the critical temperature. The departures of C_, Cy 
and the speed of sound are substantial. They are about 1% at 298.15 K, 5 - 10%, at 288.15 K 
and 20 - 40% at 283.15 K (0.8 K above critical). An analytic representation cannot describe 
the divergence of Cy and the accompanying decline to zero of the (thermodynamic) speed of 
sound while the divergence of Co is not strong enough; so these large departures are not 
surprising. The authors themselves did not claim accuracy in the critical region. Since we 
have represented the speed-of-sound data to about 1% in accuracy, we are confident that our 
model represents thermodynamic behavior, such as a weakly diverging Cys accurately. 


Ii The Effect of Pressure and Temperature Errors on Density 

In custody transfer, the density of ethylene will be calulated from its measured pres- 
sure and temperature. Since these quantities are measured with limited accuracy, an error 
in density will result that is given by 


Solem 6P 

“p = (PK) P 9; (8) 
or 

OD tS op 

5 ay 6T Ky Sit oT. (9) 


In view of the strong divergence of the compressibility Ky and the expansion coefficient On 
the density error resulting from inaccuracy in the measurement of P or T is a maximum at the 
critical point. We show the region in P-T space that should be avoided if an accuracy of 
0.1% in density is desired and the pressure is measured with an accuracy of 0.1, 0.05 and 
0.02%, respectively, (Fig. 21); or the temperature is measured to 0.02, 0.01 or 0.005 K, 
respectively, (Fig. 22). 


8. The Effect of Impurity on Density 


The present formulation of the equation of state of ethylene is based on the PYT data 
obtained for 99.993% pure ethylene, with ethane as the prime impurity. The formulation 
permits the calculation of the density of an ethylene sample from its pressure and temper- 
ature. If impurities are present in the sample of interest, the predicted density will 
differ from the actual one. It is of interest to know what effect a known impurity has on 


density, and how the effect varies with the temperature and pressure. 


Under simplifying assumptions, this question can be answered. We will assume that the 
impure fluid is in corresponding states with the pure host, and that the impure system 


8 


displays pseudo-critical constants that are shifted from those of the pure host by amounts 
linear in concentration: 


ete Tk (Ieax)) . 
Tair Re Gi bxt) es 
Meg bay N (1 + (a-b)x) (10) 


Both assumptions are reasonable and readily justifiable in dilute mixtures; the only diffi- 
culty with the model is the actual assignment of the values of a and b, parameters charac- 
teristic of the host-impurity interaction. They can be determined from the critical para- 
meters of the pure components in the case that the critical points are connected by a 
continuous critical line. For the sake of developing our estimate of the density error 
introduced by impurity, we will consider these parameters a,b as known. 

An equation for the impurity-induced density error was derived by us [3]. It was 
obtained by applying the law of corresponding states to the pure host and the dilute mixture 
at the same pressure and temperature. The following result was obtained for the change in 
density: 


[o*) 
Ko) 


-x[(a-b) + b PK. - ato, | 


+ 


o| 


= -x[(a-b) + bPK, - a(PK,) (p 2)] (11) 

where Ky is the compressibility, Op the expansion coefficient of the pure host at T, P. It 

is clear that Eq. (11) permits calculation of the impurity-induced error in density at any 
point in P-T space once the constants a, b are known. The relation (11) should not be applied 
too close to the critical point since the coefficient of the term in x diverges. More 
seriously, an expansion in x has as little justification as one in p, which is known to be 
invalid at the critical point. 

A scrutiny of Eq. (11) immediately reveals that impurity-induced density errors will be 
largest near the critical point. This is a consequence of the fact that both Ky and Op 
diverge strongly at this point. The precise nature of this density error will, of course, 
depend on the values of a and b. Since these values form the most controversial part of our 
model, we will consider only a few extreme cases. 

(1) a= -.5, b= 0: a highly volatile impurity that depresses the critical tempera- 

ture of the mixture strongly without affecting the critical pressure. 

(2) a= .5, b= 0: a highly nonvolatile impurity that increases the critical tempera- 

ture of the mixture strongly without affecting the critical pressure. 


The cases (1) and (2) have opposite effects on the density. The impurity-induced density 
error is given by 


“2 = +0.5 x [1 - (P Ky) (ey (5) 


where the + sign refers to case 1, the - sign to case 2. In Fig. 23, we show the region in P-T space 
that should be avoided in custody transfer if an accuracy in density of 0.1% is desired and if various 


amounts of an impurity very different in volatility from ethylene are present. 


oF 


Tables of Thermodynamic Properties 


In Appendix C we present the thermodynamic properties of ethylene in the critical region. It was 


not possible to organize this table along isobars, as in [5]. The extremely rapid variation of the 


"density"-like properties p, U, H, S with P and T would preclude interpolation in such a table. We 


have therefore tabulated the properties along isochores, or at constant density, using temperature as 


the second entry (Table V). In addition, we present tables of saturation properties first with tempera- 
ture, then with pressure as an entry (Table VI, VII). 


The computer program required for the generation of thermodynamic properties of ethylene is listed 


in Appendix D. Appendix E contains units and conversion factors. 
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Table I. Parameters in the Thermodynamic 
Potential for Ethylene 


Critical Exponents 


IGhined)) ia a eee B = 0.325 
6 = 4.82 
A,= 0.50 
Critical Parameters Ure = 282.3452 
©, = 7-634 mol/dm3 
Pe = 5.0403 MPa (derived) 
Parameters in 
Scaling Function aS 2S. SO) 
pie. UeSIaly 
Kk & OsoZ 
c = -0.007833 
b2 = 1.67472 
Pressure background i 
parameters Pe = 5.3439 
P, = -16.9686 
P,, = -0.181762 (derived) 
Caloric background 
parameters Uo = -36.4812 
My 2 eeV/osecy 
iy = HV. S008 
gy = 10.062 


liz 


APPENDIX A: 


A.l 


ye 


Bey! 


et 


Cc 
Vv 


Revised and Extended Scaling Equations for the Thermodynamic 
Properties of Fluids. 


Reduced thermodynamic quantities 


Le Fee, cre er Le 
The aie iT; P : RTP : 
c c 
oe ja Oe ef isveane 
d 3 a 
P. i VA Va Be 
(A.1) 
A Sue fee ac ait ee 
MP Vane oP? 2S OnE |i 
¢ € U fr 
Ce LT ec: Tt 
mV aC (Cu ee tee 
Vie Pe, P av SP 
c Cc 


(T is temperature, wp is chemical potential, P is pressure, p is density, 


U is energy, S is entropy, A is helmholtz free energy, H is enthalpy, 


V is volume, Cy is heat capacity at constant V, CS is heat capacity at 


constant P). 


A.2 


S 


E 


Thermodynamic relations 


UdT + pdy 
-UdT + jdo 
-TdU + pdu 


-TdU - udp 
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(A.2) 


A.3 Fundamental equations 


Au =u - (7) 

P = Pi) + Ap + P, AAT + AP 
with 

s AS oan ~ 5 

He CD) eee ete en CATy) 


A.4 Derived thermodynamic quantities 


ce 


1 + P AT + OAP 
el dAL/ .~ 
AT 
des canis Gna 
Siete Ep Au + = 
dT dT 1] OAT ne 
u 
97 Ap 
ND = 
oA / aT 


(A.3) 


(A. 4a) 


(A. 4b) 


(A.5) 


(A. 6a) 


(A. 6b) 


(A.7) 


(A.8) 


(A.9) 


A.5 


with 


A.6 


eee As [yee 2.1, /aAP UNS 

iar) a da Dg oN age  -¥ DALaAT 
6 Xp Ai 40 Xp P4ueAT 

BN. ro), con festa 

7? aT? : ale 6 a 


Critical exponents 


ec | Rye 8 oA, 
ee yer 2 OS 
ct= e(6rl)) | Skah ¥ = 6 (6-1) 


Parametric equations for singular terms 


in a -Pan ole 


AT = r(1-b76*) - cAji 


(A.10) 


(A.11) 


(A.12) 


(A.13) 


(A.14) 


(A.15) 


(A.16) 


Le ple eee 
AP = 5 r* “hak p. (6) (A.17) 


(35) = ; [=F2x;0 + cr 814k, 55 (6)] (A.18) 

Gal z 2 ri 4ak,5, (8) (A.19) 

(o) : 2 fe u, (8) +2er"i"Tk,v, (8) 1or Mak, w, (0 (A.20) 
dA he 1=0 Fn 11 ivi 

eu - Ey [8 hey, (0 + cx ak, (0) (A.21) 

(23) : Es r ak,w, (8) (A. 22) 


A.7 Auxiliary functions 


P3(8) = Py; + P38 + Py38 (A.23) 
s, (8) = s,, + Se : SU) = 2a a (A.24) 
HO) Go © EGR eSes 6 BORA e” (A. 25) 
u, (0) = [1 - b”(1-28,)07]/q(8) (A. 26) 
v, (0) = [B,(1-367)@ - B8(1-67)6]/4(8) (A. 27) 
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w, (0) = [(1-a,)(1-307)s, (8) - 86(1-0°) 6s! (8)}/q(8) (A. 28) 


with 


2 
B6-3B;-b a.Y; 


Fog maar are 
2b (2-a;) (1-a, Ja; 
Z 
B6-3B.-b a, (286-1) 
Da ae ce (A. 29) 
2b (l-a. )a. 
ivi 
if 286-3 
Pai 2a. 
i 
Soi * (2-05) Poi 
65-38 
I eS 
2i 2b70, (A. 30) 


A.8 Two-phase properties 


Variables: 
(6) = geil 
Au = 0 
~ 2 
AT = r(1 - b) (A.31) 


Vapor pressure: 


1 
Pap = Pi (7) + ae r i ak. p, (1) (A.32) 


WW 


coexisting densities: 


pep 1 
L+' Vv = = 4 -ay 
ae 1 + Pay AT + ito (S93 ak;s; (1) 
p,-#. 1 
Lavine Bi 
> io rk; (A.33) 
helmholtz free energy: 
} 2 
Sis Sieve yc ~ “A. | 
APSO UNC) IPE AC) as ee ro i 'ak;p,; (1) (A. 34) 
energy: 
dP lain) 1 Me 
U= capes i Some 0) »> (2-@.)r 1 ak =P; (1) (A.35) 
dT dT l=by 250 
entropy: 


1 = 1 
Ay ei A ol) os ne ap) (A388) 
i=0 oat Qin rai i at ot 
1-b 1i=0 
specific heat Cc, 

| athe “herdon Ne “a, 

Spe epi Ok eon | a (2-a,) (1-a,)r ak-p,; (1) (A. 37) 

Tr oat dT 1-b°/ i=0 7 
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Appendix B. Comparison with Experiment 


In Table II, the experimental PVT data are compared with the formulation (a) as pressure differences 
under the assumption that the densities are exact and (b) as density differences under the assumption 
that the pressures are exact. In Table III, saturation pressures and densities from two sources are 
compared with the formulation. In Table IV, the speed-of-sound data [4] are compared with the 
theoretical prediction. The density (DENS) is in mol/dm? and the speed of sound w in m/s. The legend 
for the various data sources is as follows 


1 DH [10] 
2A HLSB [3], 99.993% sample 
2B HLSB [3], 99.999% sample 
3. TwW [11], 
a IEEE 
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kabler lle 


TEMPERATURE 


298. 
293%: 
. 1500 
.6500 
.3440 
.1393 
.6394 
.3394 
.0354 
- 7394 
3.1500 
.6500 
3.6500 
-3440 
gp USI} 
0.8393 
-7393 
. 1509 
. 1500 
&.1500 
2.6500 
.6500 
. 3440 
.8393 
. 7893 
. 1500 
. 1560 
9.1500 
-6€500 
. 6500 
. 3440 
. 1893 
.0993 
.0949 
.0843 
.0936 
.0893 
.0891 
.0877 
. 1500 
-6500 
.3440 
.2193 
.1916 
.1911 
. 1909 
1893 
.17%8 


K 


1500 
1500 


Comparison of the experimental PVT data 


from four sources with the formulation 


GAMKIninnngyqynqgnggavwiwgagnwnnaiwgqgnanrtAaAhHUUUOHONIHAAAMUGTNIAAAAUAIONA 


PRESSURE 
EXP 


32377 
33162 
.54857 
.49438 
.24955 
.02650 
.93481 
-88010 
.82552 
»77146 
.89748 
.68640 
. 34587 
£12725 
-92924 
. 88062 
.86455 
-46225 
.68170 
.O06666 
.53899 
» 24431 
.05644 
.98547 
»97849 
.29472 
-56545 
.99031 
-49642 
. 22054 
.04511 
.02483 
-01308 
.01251 
.01246 
.01233 
.01181 
.01176 
-01156 
.23391 
. 21388 
-04288 
.02695 
. 02346 
-02344 
02335 
-02316 
.C2097 


(MPA) 


VIIa nrIaninainngIngawagnogoansAAbhOUNOO WTA HAHAH UNHADAALNOINIUON DD 


CALC 


. 31968 
. 32880 
. 54693 
-49414 
. 24957 
.02643 
-93461 
.87976 
. 82505 
.77051 
89935 
. 68744 
- 34697 
12824 
.92977 
. 88086 
.86459 
-46412 
. 68233 
.06664 
. 53879 
24431 
. 05660 
.98552 
.97853 
. 29582 
56539 
93965 
- 49573 
. 22022 
-04512 
.02482 
.01306 
-01249 
.01241 
.01232 
.01176 
.01173 
-01155 
. 23474 
Honisioe 
.04269 
. 02697 
.C2345 
02338 
0233: 
O88 1S 
.0208) 
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DENSITY 
EXP 


10. 
NOR 
10). 
VO 


MOoODMDMEOMDOBDMDOMAWDAOAMDOAMDAAAAMAMNINNONONWNNHOHOWOL 


542100 
544500 
546700 
549500 


VBS 92010 
~556590 
.551900 
.551200 
.551300 
ADS) 1) SONG) 
025960 
.029200 
.030300 
-031500 
AOSZNCO 
.032200 
. 032300 
. 285400 
. 287600 
. 289400 
. 290000 
. 291300 
. 292400 
. 292500 
. 292600 
.910500 
. 912600 
.914300 
.915800 
.916700 
.917200 
.917300 
.917300 
.917300 
.917300 
.917300 
.917300 
.917300 
.917300 
. 757700 
. 763800 
. 764300 
. 764300 
764400 


764400 


. 764409 
. 76440" 
. 764404 


(MOL 


DODDDDDDDMODNDDDDHDMDMHDMDTMDDOWKVOKVHOHO 


/DM3 ) 
CALC 


547447 
548586 
.549579 
550146 
550143 
551053 
551397 
.552670 
553901 
.556070 
021107 
026154 
025659 
.025416 
.027673 
029998 
.0319c0 
. 281344 
285694 
289482 
292384 
.291775 
. 288973 
. 230874 
291361 
.907691 
912840 
918022 
923084 
923166 
.S16409 
.917943 
.918628 
918978 
.921296 
.918123 
921527 
919558 
.918114 
755431 
772971 
763615 
. 762402 
766152 
771465 
. 753508 
705213 
784795 


TEMPERATURE 


303. 
298. 
293; 
289. 


K 


1500 
1500 
1500 
1500 


.6500 
.6500 
.6500 
. 1500 
.6500 
.5000 
. 3440 
.1500 
. 1500 
. 1500 
. 1500 
.6500 
4.6500 
.6500 
.1500 
.6500 
. 3440 
. 1500 
.1500 
.1500 
. 1500 
.6500 
.6500 
.6500 
. 1500 
.6500 
.3440 
. 1500 
511 SOKO) 
.1500 
. 1500 
.6500 
.6500 
.6500 
. 1500 
.6500 
2.3440 
. 1500 
. 1500 
.1500 
. 6500 
1E509 
SG 
. 3440 


0392 


“J 


LHOOYUKVUVUOPOOHUWIAGHHTIWIWIH OV 0 VW U1 OO OD O SN OF OF OF OF OF UT U1) M ~I OF UF OF UT OF UO MI GO) @ 


Table II. (Continued) 


PRESSURE 
EXP 


-47281 
.88257 
.29461 
. 82681 
- 42007 
- 30464 
- 18952 
.13214 
-07498 
-05789 
-04020 
- 45336 
. 86844 
. 28553 
-82149 
-41783 
- 30310 
18875 
.13169 
-07481 
.04007 
- 43832 
.85752 
-27857 
.81749 
-41617 
. 30198 
- 18817 
13139 
-07473 
.04016 
-40673 
-83445 
. 26383 
. 80890 
»41252 
- 29966 
. 18698 
.13070 
.07449 
.04011 
.58548 
-10120 
ot it 7s} 
363233 
.16997 
. 38762 
.02543 
- 92780 


(MPA } 


HLYDAMHUMDOUNHUTUTMHUUNAODYNDATITUOUNTHUTVADDVGMTUTHOKTUHOOVNNMNMTHUTUIADN 


CALC 


- 47347 
. 88269 
- 29437 
. 82639 
- 41985 
. 30441 
. 18940 
.13210 
.07498 
.05789 
.04016 
- 45393 
. 86851 
. 28531 
VO2Nt7, 
.41766 
. 30300 
. 18869 
suse 
.07486 
.04016 
- 43893 
.85760 
. 27834 
.81716 
.41598 
. 30191 
.18815 
.13141 
.07477 
.04016 
-40710 
. 83453 
. 26359 
. 80863 
- 41240 
. 29959 
. 18699 
. 13076 
.07457 
.04016 
.58491 
. 10096 
ef WAS) 
. 36860 
. 17027 
. 38666 
-02552 
199739) 
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% 
DIFF 


009 
.002 
004 
.007 
004 
.004 
.002 
.001 
.000 
.000 
001 
.008 
.001 
.003 
005 
.003 
.002 
001 
.000 
.001 
.002 
.008 
.001 
004 
.006 
.004 
.001 
.000 
.000 
.001 
.000 
-005 
.001 
-004 
.005 
.002 
.001 
.000 
.001 
.002 
001 
.009 
.004 
00} 
.005 
.006 
.015 
.002 
.004 


GOTOH DAMOVIVAVVIAINVYVIINIVYVNYV YIN YIN YI YIN YI YIN IS NINN NNN NINN 


DENSITY (MOL/DM3) 


EXP 


.672200 
.674000 
.675800 
.677200 
.678500 
.678900 
.679200 
.679400 
.679600 
.679600 
.679700 
.621770 
.623570 
.625370 
.626800 
.628060 
. 628420 
.628780 
. 628950 
.629130 
.629240 
. 582800 
. 584500 
. 586300 
- 587800 
5989000 
.589400 
.589700 
.539900 
.590100 
,590200 
- 499400 
.501200 
. 502900 
- 504200 
. 505600 
. 505900 
. 506300 
. 506500 
. 506600 
. 506700 
.579700 
.531300 
. 582500 
. 583600 
584200 
899000 
. 903300 
. 903460 


GUONTODOBOMOVWYIVVINVIYVSINYINS YG YIN VINNY YIN NY SIV YIN NIN YIN NNN NN 


CALC 


.670365 
J6735E°S 
.677184 
- 681333 
.683527 
.687088 
.687745 
.684553 
.67E838 
.675203 
. 762056 
.620178 
.623313 
.626623 
.629959 
.631947 
.6321914 
.632357 
.626546 
. 606356 
-492949 
.581088 
.584204 
LOS: 
.591089 
.593476 
.591987 
.591210 
.587625 
Noa 2O 2 
. 534842 
.498349 
.500903 
. 504331 
,505977 
.508351 
. 508313 
.505527 
- 428144 
.471538 
. 392577 
.581888 
.982547 
. 582234 
.578960 
.573036 
.902252 
. 900275 
.897413 


me ek es a as ot es 


TEMPERATURE 


298 


260. 
209% 


279 


298. 


282 
282 


282. 


282 
282 
282 
282 
282 


282. 


282 
282 


282. 


282 


282. 


ik 


.1500 
0394 
5394 
. 3625 
1500 
.3700 
.3440 
2143 
. 1893 
. 1500 
.1500 
.1500 
.1500 
.6500 


.6500 
. 6500 


. 1500 


.6500 
.5500 
-4500 
. 3440 
. 1500 


. 1500 


. 1500 
-6500 
- 6500 
. 2446 
. 1500 
. 1500 
}, 1500 


.6500 


.6590 
2.5419 
2393 
. 1393 
9393 
.8393 
7393 


.3500 
3600 
.3700 
3800 
$3900 
2000 
4100 
4200 
. 4300 
4500 
5500 


Table II. (Continued) 


PRESSURE 
EXP 


39638 
-84374 
a ehtAS) 
- 71942 
- 23094 
-04648 
.04315 
.02663 
.02347 
~44551 
. 86296 
28229 
.81975 
-41734 
. 30307 
. 18886 
.13189 
.07514 
.06378 
.05243 
-04045 
- 56656 
-08870 
- 70401 
. 36469 
. 16848 
-03822 
. 38612 
-96311 
.62133 
. 31829 
~14244 
.02570 
.01624 
-00720 
-98901 
-97992 
.97071 


VASAT ATUOUMTIITIUOVOOMUTUNUTUTATUTITOGANUTUGTUN Hh DAA 


HAD 


5.04197 
5.04288 
5.04403 
5.0451€ 
5.04626 
5.04744 
5.04858 
5.04978 
5.05083 
5.05324 
5.06452 


(MPA) 


DHSBANYVUNNUTUIOWUTUUVHDOUHUUKTNKHHATIHIOODNUNUHUNH HHO 


U1 OT OF OF OT OT OT OF 1 UT 


CALC 


. 38978 
- 84231 
- 74989 
- 71724 
- 23176 
04613 
.04280 
.02627 
02309 
. 44643 
. 86306 
. 28184 
.81916 
- 41682 
. 30245 
. 18842 
. 13156 
.07482 
06349 
.05216 
.04016 
. 56593 
.08818 
. 70374 
. 36453 
. 16836 
.03807 
TOBSny, 
- 96213 
- 62040 
. 31759 
. 14189 
.02534 
.0159% 
. OC6E8 
. 98878 
.97979 
. 97062 


-04084 
.04198 
.04311 
.04424 
.04538 
.94651 
.04764 
.04878 
.04991 
.05218 
06352 
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.018 
021 
.020 


DENSITY (MOL/DM3) 
EXP 


10. 
10. 
10. 
10. 


GDC U UU UT UU ODNDONDONVYANNIYINYN YY YNNOOWAOM 


NNYVVN SIN NNN NS 


617200 
626500 
626800 
626800 


. 750000 
. 756600 
. 756600 
. 756600 
. 756700 
. 602300 
.604100 
.605900 
.607300 
.608500 
.608900 
.609300 
. 609400 
. 609600 
.609600 
.609700 
.609700 
.511200 
.51280C 
.51400C 
.515000 
.515600 
.516000 
.894200 
.895500 
. 896600 
.8976500 
.898100 
. 898500 
.898500 
. 898600 
.898600 
.8938600 
.898700 


.635000 
.635000 
635000 
.635000 
.635000 
635000 
635000 
.635000 
.635000 
.635000 
.63500C 


CALC 


10.624430 
10.632286 
10.635088 
10.636740 
.747744 
.784006 
. 784956 
.794163 
. 797239 
.599740 
603709 
.608538 
613134 
620581 
.631131 
640665 
.652152 
. 739052 
787493 
.848875 
. 732198 
.513592 
515674 
.516311 
.517597 
.519788 
534596 
.897429 
.900023 
.902760 
.905107 
. 907296 
. 908266 
.908148 
. 908367 
.906675 
54906706 
5 ,902345 


ATHAKUMAIKTUHIUADNDINIDMDMOHDVYVVNVNANN YANN NNN OOODWD 


. 797934 
. 800451 
809395 
.818024 
.827379 
.836230 
.845222 
.854504 
. 862487 
.878148 
. 930706 


NYVNNNNNSNN ANN 


EXPT 
a) 


aauwewnrnnenn neo wt tre were er ew ewer w tw rt tr wrrwrorstwrerer ror stsrststst—s 


TEMPERATURE 


282. 
282. 
282. 
282. 
283. 
283. 
283. 
284. 
285. 
286. 
XH 
288. 
293. 
298. 
303. 
283. 
288. 
293. 
298. 
303. 
283. 
238. 
293. 
298. 
303. 
283. 
288. 
293. 
298. 
303. 
283. 
288. 
293. 
298. 
303. 
233. 
283. 
293. 
298. 
303. 
283 
288. 
293. 
298. 
303. 
283. 
288. 
293. 
298. 
303. 
283. 
288. 
293. 


K 


6500 
7500 
8500 
9500 
0500 
1500 
3500 
1500 
1500 
1500 
1500 
1500 
1500 
1500 
1500 
1500 
1500 
1500 
1500 
1500 
1500 
1500 
1500 
1500 
1500 
1500 
1500 
1500 
1500 
1500 
1500 
1500 
1500 
1500 
1500 
1500 
1500 
1500 
1509 
1500 
1500 
1500 
1500 
1500 
1500 
1500 
1500 
1500 
1500 
1500 
1500 
1500 
1500 


OUNUDAINOMKVIIVIAVAUVTHIVOOTIVTIION HGH HTIMDHSOTVUIMNAMUTTGOOUMTUTHIOMHMTHUNAAU HO OM WI OI 


PRESSURE 


EXP 


.07586 
.08723 
.09858 
- 10994 
12136 
13277 
. 15567 
. 24722 
. 36197 
47701 
.59247 
. 70808 
. 28886 
-87365 
-46149 
. 10345 
-54964 
- 98672 
- 42028 
-84890 
.11928 
. 60692 
.08774 
»56526 
.03963 


12707 


.65366 
.17782 
. 70124 
22319 
-13169 
.69652 
. 26526 
-83652 
-40872 
13594 
. 74088 
.35743 
-98029 
-€0625 
-14271 
- 79342 
- 46250 
.14174 
.82674 
»15732 
.86313 
-59227 
- 33468 
.08490 
- 18994 
-96313 
. 76196 


Table II. (Continued) 


DOMNUDIBDUTIHIVODVTIUVIAHDIGIIVDOVTKVAIDSUMTIUIVDMNDUVUHNDWWNWNY DO VW OF VW VI OF UI OF OF OF OF OT 


(MPA) 


CALC 


.07488 
. 08624 
.09761 
. 108983 
- 12037 
UCU LAs) 
. 15455 
24593 
. 36053 
-47548 
.59074 
- 70626 
. 28704 
.87171 
. 45904 
- 10241 
.54816 
.98515 
-41773 
. 84738 
.11846 
. 60572 
.08577 
.56281 
.03798 
.12625 
. 65253 
.17599 
. 69893 
. 22170 
TSO 
. 69483 
. 26308 
83420 
. 40733 
. 13481 
. 73870 
. 35480 
. 97763 
- 60492 
~149141 
- 79106 
- 46020 
a Vee 
. 82622 
. 15588 
. 86107 
.59056 
. 33347 
.08531 
. 18874 
.96149 
. 76038 


28) 


9 
% 


DIFF 


.019 
019 
019 
.019 
.019 
020 
.022 
.025 
.027 
.028 
.031 
.032 
.029 
-028 
033 
020 
.027 
.026 
.040 
.022 
.016 
021 
032 
.037 
.023 
.016 
-020 


030 
-034 


021 


.019 
.030 


.035 


034 
JOS 
.022 
938 
.041 
.C38 
Oth? 
.025 
041 
036 


.028 
-007 
.028 
.035 
.026 
SOU C 
.005 
.023 
.028 
O23 


OWWNDNHNHOWHUDAWMDWDAMDADAMWMDWNAIAIANYVAINNINNVONDNDDNDADDADADADADYNNYVYNNNYNNNNYNN SN 


DENSITY (MOL/DM3 ) 


EXP 


. 635000 


. 635000 


635000 
. 635000 
.635000 
635000 
. 635000 


- 6355900 


.635000 
.635000 
.635000 
. 635000 
.635000 
.6350090 
. 635000 


. 000COO0 


. 000000 
. 000000 
. 000000 
. 000000 


. 500000 
. 900000 
. 500000 
. 500000 
. 500000 
. 000000 


. 0O00CO0O 


. 000000 


. 000000 


000000 
.500C00 
.500000 
500990 
500000 
. 500000 
. 000000 
000000 
.0C0000 
. 000000 
.000C00 
. 500000 
500000 
. 500000 
.500000 


. 500000 
. 000000 
.000000 


. 000000 
. 000000 


. 000000 
. 500000 
. 500900 
. 500000 


OWWUWDWOWHWUOMDAWDMDAWDAAADAYININNYNAININININNNANMNVDADAADONNNNNYANNNNNNNNNN 


CALC 


.914154 
873055 
. 830066 
.795675 
.776276 
759454 
741249 
.696503 
.675567 
.664942 
.660746 
.656658 
645534 
.642283 
.641778 
023921 
.011648 
.007569 
.008853 
.004152 
539665 
.512006 
.510856 
.509280 
.504772 
.080309 
013259 
.010858 
.009055 
.C04380 
624604 
.52041C 
.512838 


. 508388 
. 503930 


. 103823 


.0235C6 
.014007 
.609298 
.002442 
.562742 
.519981 
.510251 


. 505906 
.501175 
.031780 
.012601 
.006005 
.0023977 
.999225 
.512497 
. 506879 


.504175 


DIFF 


.656 


.850 
.630 
392 


.392 
Sz 
,2a4 
.138 
.095 
.089 
.399 
.194 
£126 
.148 
.069 
.610 
.185 
.167 
.143 
.073 
NAAT. 
.189 
.155 
.129 
.063 
.501 
a 272 
aT 
.119 
.052 
. 296 
.294 
2175 
SING 
.043 
.738 
.235 
a2 
.069 
.014 
Fo5D3 
.140 
.067 
033 
.009 
32 
.072 
.044 


EXPT 


~ 
Oo 


Sa a a CC Sl ee ee ee 


TEMPERATURE 


298 
303 
223 
288 
293 
298 
303 
283 
288 
253 
2S8 
303 
280 
280 
260 
280 
260 
286 
280 
280 


K 


. 1500 
. 1500 
. 1500 
. 1500 
. 1500 
. 1500 
- 1500 
. 1500 
. 1500 
. 1500 
. 1500 
. 1500 
. 1500 
. 1500 
. 1500 
1500 
1500 
. 1500 
. 1500 
.1500 
. 1500 
. 1500 
. 1500 
. 1500 
. 1560 
.1500 
. 1500 
- 1500 
.1500 
.1500 


. 1500 


. 1500 
. 1500 
. 1500 
.6500 
.6500 
-6500 


.6500 
-6500 
-6500 
.6500 
-6500 
.6500 


.6500 
. 1500 
. 1500 
. 1500 
. 1500 
. 1500 
. 1500 
. 1500 
. 1500 
. 1500 
.1500 


VUOUUCTIUTUUUO BS HPHBAHBHHHADHDSDSHADAHDHAHHAAHAUHBDAAAHDAHAODNOVNAOVNPVOUAN 


Table II. (Continued) 


PRESSURE 


EXP 


-57508 
- 39804 
.25792 
- 11098 
.99048 
.88281 
- 79026 
. 38690 
33296 
. 30284 
. 29048 
. 28873 
-83757 
.81814 
-81214 
. 80596 
. 80263 
.79995 
. 79989 
.79270 
- 79960 
.01853 
-92924 
.91885 
.91325 
-91168 
.91022 
- 90938 
. 90885 
- 90837 
- 90826 
- 90822 
. 90804 


° 


90601 


-97260 
-96918 
. 96762 
. 96606 
-96472 
. 96393 
- 96300 
- 96067 


95740 


.94814 
.20165 
.09211 
.04078 
02694 
02425 
02198 
.02061 
.01990 
.01959 
01951 


OQOnIannain»niwnsaHLHPHDHDHHDDHAHAHAHHAHADAHAHAYHDHAALHBHDALHSANOWDNAIADNMNANIAONWAN 


(MPA) 


CALC 


57459 
.39918 
25693 
10974 
.98814 
.88257 
. 78816 
. 38340 
32812 
. 29650 
. 28063 
27603 
83427 
.81532 
. 80947 
80366 
.80031 
.79910 
.79910 
.79910 
79910 
01542 
.92778 
.91745 
.91176 
.91022 
.90873 
.90777 
90731 
. 90736 
.90736 
. 90736 
. 90736 
-90524 
.97138 
.96794 
. 96647 
.96479 
96353 
. 96262 
. 96253 
95983 
. 95639 
.94693 
19852 
09120 
.03979 
02584 
02305 
.02096 
.01945 
01884 
01845 
.01843 


24 


ves) 


—_ «A 
DOMDDHHWHOOUUWUMTHOWUWUDNHNWHOUUHDDNONLeCUOUUODVNOWUODYNDAYH 


DENSITY 
EXP 


.50C000 
. 500000 
.000700 
.0009000 
.900009000 
,000000 
.000000 
.500900 
.5006C0 
. 500000 
. 500000 
. 500000 
. 500900 
- 400900 
. 365500 
. 329100 
.397900 
. 900000 
. 250900 
. 8009090 
.850000 
. 500000 
. 000500 
.899000 
- 833900 
.813400 
. 793600 
. 780300 
. 773700 
. 767300 
. 300000 
. 100000 
. 200000 
. 500000 
.592500 
. 530000 
. 500000 
.462700 
. 432000 
. 408050 
. 900000 
. 800000 
. 700090 
. 500000 
.500000 
.000000 
.500C00 
. 200000 
. 100000 
.000000 
.950000 
.845200 
.802400 
. 800C00 


(MOL/DM3 ) 


_a 


HODDOWONONDOOCOUUNUNUNOMOOHONNDHOOUHADWOOCHHONHNNDTOUDe! 


a 


CALC 


. 500559 
.498235 
.005284 
003532 
.004588 
.000366 
.002596 
.510550 
.509354 
.509159 
.511444 
.512414 
.515588 
.415616 
381149 
.343470 
. 322035 
.905268 
253402 
.801760 
.801668 
.512315 
.012537 
.913674 
.850774 
.831534 
.612843 
.601638 
.794532 
.781115 
.311185 
.108141 
. 203784 
.51€828 
.612074 
.552898 
.522903 
.490067 
.459807 
.440300 
.904845 
. 829248 
.727739 
521962 
.510703 
. 006037 
.515163 
.231517 


143214 


048220 
968625 
.918681 
.888712 
883761 


EXPT 


—e ee ek me me sk we wk ek we a ek kk tt wk tt wt wk ww wt kw ath wt wt tt wt wt se es wt ss wt 


TEMPERATURE 


282. 
282. 
282. 
282. 
282. 
282. 
282. 
282. 
282. 
282. 
232. 
282. 
282. 
262. 
. 2500 
282. 
.2500 
. 3500 
. 3500 
. 3500 
282. 
282. 
282. 
282. 
282. 
282. 
2&2. 
282. 
. 3509 
282. 
282. 
282. 
282. 
- 4500 
4500 
.4500 
2282. 
- 4500 
.4500 
262. 
282. 
282. 
262. 
282. 
282. 
282. 
.€500 
6500: 


282 


282 
282 
282 
282 


262 


2382 
282 
282 
232 
282 


282 


282. 
282. 
282. 


282 
282 


K 


1500 
1500 
1500 
1500 
1500 
1500 
1500 
2500 
2500 
2500 
2500 
2500 
2500 
2500 


2500 


3500 
3500 
300 
35090 
3500 
3500 
3500 
$500 


3500 
£000 
4000 
4500 


4500 


5500 
5500 
€590 
6500 
6500 
6500 
6500 


6500 
6500 


- 6500 
-6500 


MM OFT OF OT OF OT OF OF OT OF OF LT UT OF UT OF OF OF OF OF UT OF OF OF OF UF OF OFT OF OF OO OT OF OF & OT OF 0 0 OF OF OFT OF B I UI I OT 1 


Table II. 
PRESSURE (MPA) 
EXP CALC 

.01843 5.01759 
.01741 5.016568 
.01606 5.01515 
01181 5.01089 
00892 5.00790 
.00530 5.00423 
.98973 4.98840 
.03540 5.03427 
.03228 5.03123 
.03142 5.03035 
.03089 5.02978 
.03035 5.02954 
02995 5.02904 
.02939 5.02843 
.02842 5.02755 
.02118 5.02011 
.99803 4.99666 
. 23803 5.23536 
.12506 5.12417 
.07067 5.06935 
.04902 5.04762 
04295 5.04166 
.04200 5.04087 
.04200 5.04085 
04194 5.04084 
04186 5.04082 
.04178 5.04067 
03°66 5.93653 
.02051 5.02931 
60613 5.00491 
04746 5.04651 
04741 5.04650 
.05343 5.05247 
.05322 5.05225 
.95309 5.05219 
.05309 5.05216 
.05303 5.05209 
.05283 5.05183 
.05235 5.05146 
.06449 5.06354 
.06446 5.06347 
.08928 5.08795 
.07993 5.07874 
.07695 5.07585 
.07616 5.07515 
.07590 5.07491 
07575 5.97479 
07556 5.07456 
.07481 5.07389 
.07383 5.07294 
.07221 5.07136 
06950 5.06867 


25 


(Continued) 


% 


DIFF 


HOM? 
LOU 
.018 
.018 
.620 
.021 
.027 
O22 
.021 
021 
022 
016 
.018 
.019 
017 
.021 
SOL 
051 
= OMthz, 
026 
.028 
.026 
022 
923 
O22 
.021 
.022 
992 
.024 
024 
019 
.018 
.019 
019 
.018 
.018 
.019 
.019 
.018 
.019 
.019 
.026 
.023 
O22 
.020 
.019 
019 


0205 


.018 
.018 
-017 
.016 


—_ 


DONVWAINNNVNDAAON I YuUNNININVOINTIANANDININNDDOWOVOVOWOUANAAAANAAMAAWHDNOUUUANANAD 


DENSITY 


EXP 


. 400000 
. 3090000 
. 200000 
.000000 
. 900000 
. 800000 
. 500000 
.O0CO000 
. 800000 
. 700000 
. 600600 
. 750000 
.600000 
.500000 
. 400000 
. 000000 
.5900000 
.500000 
.000000 
. 500000 
. 000000 
.500000 
.000000 
. 750900 
. 500000 
. 250000 
.OCODCO0 
.500000 
.000C909 
.500006 
.€50000 
.600000 
. 000000 
. 750000 
.650000 
. 600000 
.500000 
.250000 
. 000000 
-650500 
. 6CO000 
. 000000 
- 500000 
.000000 
. 750000 
. 650000 
.600000 


. 500000 
. 250000 
.000000 
. 750000 
. 500000 


(MOL/DM3 ) % 
CALC DIFF 
6.454023 - .844 
6.390674 -1.439 
6.265869 -1.062 
6.036959 -.616 
5.932741 - 555 
5.827858 - .480 
5.520626 -.375 
9.047636 -,529 
8.872435 - .823 
8.792945 -1.068 
8.716073 -1.350 
6.81894S -1.021 
6.701767 -1.542 
6.785221 -4.388 
6.515809 -1.81C 
6.040185 -.67C 
5.520597 -.374 
10.508902 - 085 
10.005617 - .056 
9.518329 -.193 
9.052224 - .580 
8.635545 -1.595 
8.176168 -2.202 
7.918312 -2.172 
7.€56728 -2.090 
7.397917 -2.0406 
7.158214 -2.260 
6.704560 -3.147 
6.042159 -.703 
5.517730 - 322 
7.851210 -2.630 
7.801216 -2.648 
6.202698 -2.534 
7.982894 -3.005 
7.892216 -3.179 
7.848177 -3.265 
7.762555 -3.502 
7.579103 -4.539 
7.925805 13.226 
7.938433 -3.770 
7.906237 -4.029 
9.039806 e442 
8.593304 -1.098 
8.191538 -2.394 
7.997926 -3.199 
7.925024 -3.595 
7.886525 73.771 
7.827609 -4.368 
7.625743 -5.183 
7.275154 -3.931 
6.882392 -1.961 
6.567970 -1.046 


Table II. (Continued) 


EXPT TEMPERATURE PRESSURE (MPA) % DENSITY (MOL/DM3 ) % 

ID K EXP- CALC DIFF EXP CALC CIFF 

1 282.6500 5.05791 5.05683 .021 6.0C0000 6.031517 > gg) 
1 282.8500 5.10045 5.09937 .021 8.000000 8.145839 - 1.823 
1 282.8500 5.09920 5.09812 021 7.750000 7.941939 -2.477 
1 282.8500 5.09867 5.09767 020 7.650000 7.846700 OA) 5 SY |} 
1 282.8500 5.09845 5.09746 .019 7.600000 7.803523 -2.678 
1 282.8500 5.09804 5.09702 .020 7.500000 7.720402 -2.939 
1 282.8500 5.09675 5.09584 .018 7.250000 7.447809 -2.728 
1 282.8500 5.09507 5.09431 .015 7.000000 7.122422 -1.749 
1 283.1500 5.13368 5.13266 -020 7.750000 7.867914 of] 827 
1 233.1500 5.12950 5.12868 .016 7.250000 7.348742 -1.362 
1 283.1500 5.12379 5.12302 SONS) 6.750000 6.804352 - .805 
1 283.3500 5.18470 5.18327 .028 9.000000 9.028782 - .320 
1 283.3500 5.16820 5.16669 .029 8.500900 &@.561573 -.724 
1 283.3500 5.15981 5.15853 025 8.000000 8.094200 = tis 
1 233.3500 5.15692 Sets 573 023 7.750000 7.856452 -1.374 
1 283.3500 5.15419 5.15318 .020 7.500000 7.598528 -1.314 
1 283.3500 5.15144 5.15055 LOMA 7.250000 7.332635 -7.140 
1 283.3500 5.14829 5.14749 .016 7.000000 7.061968 - 885 
1 283.3500 5.14440 5.14357 .016 6.750000 6.798297 -.716 
1 283.3500 5.13908 5.13827 .016 6.500000 6.533322 - .520 
1 283.3500 5.12159 5.12056 .020 6.0% 0000 6.021887 - .365 
1 283.2500 5.08842 5.08697 .028 5.500000 5.516559 - .302 
1 284.1500 5.57352 5.56985 . 066 10.500000 10.509930 - .095 
1 284.1500 5.42568 5.42462 .020 10.0C0000 10.004796 - ,.048 
1 284.1500 5.34153 5.34013 .026 9.500CC:0 9.511175 -.118 
1 284.1500 5.29548 5.29375 .033 9.000000 9.025101 -.279 
1 284.1500 5.27042 5.26868 .033 8.500000 8.544784 D5 BAU 
1 264.1500 5.25551 5.25395 .030 8.000000 8.062361 -.780 
1 284.1590 5.24964 5.24834 .025 7.750000 7.809672 -.770 
1 284.1500 5.24427 5.24314 .022 7.5G0000 7.554759 -.730 
1 284.1500 §.23885 5.23789 .018 7.250000 ’. 294869 -.619 
1 284.1509 5.23306 5.23214 .017 7.000600 7.038093 -.544 
i 284,150 5.22627 5.22541 .017 6.750000 €.779889 - 443 
1 264.1500 5.21797 5.21709 O17 6.500000 6.524163 -.372 
1 234.7500 5.19400 5.19278 .023 6.000060 6.020626 - .334 
1 284.1500 5.15364 5.15210 030 5.500000 5.515228 25400 
1 5€ 1500 5.57503 5.67387 -020 7.250000 7.264667 -.194 
1 285.1500 5.63106 5.63007 .018 6.720006 6.760909 -.162 
4 283.1700 5.10770 5.10684 .017 6.060000 6.081330 o 5 She) 
4 287.3700 5.48740 5.48530 .020 6.060006 6.069942 -.164 
4 293.3400 6.01530 6.01471 -010 6.0€0G00 6.062838 - .047 
4 298.1800 6.43890 6.43815 .012 5.056000 6.058623 - .043 
4 285.5400 5.42540 5.42420 -022 8.063000 8.037968 -.310 
4 291.9000 6.21100 6.21070 .005 8.063000 8.064815 -.023 
4 292.7500 6.31690 6.31702 .002 8.063000 8.062361 .008 
4 293.1500 7.20000 7.20188 .026 10.362000 19.359027 .029 
4 293.6900 7.30330 7.30472 .019 10.362000 10.359812 021 
4 298.1400 7.99060 8.C1997 . 368 10.190000 10.148722 -405 
4 298.1400 6.99690 6.99828 .026 8.072000 8.067244 .059 
4 298.1400 6.69610 6.69977 .055 7.007000 6.992544 . 206 
3 283.6500 5.14322 5.14226 .619 5.904244 5.920517 >. 276 
3 285.6500 5.31888 5.31822 .012 5.904244 5.911351 -.120 
3 289.1500 5.62216 5.62158 .010 5.904244 5.998078 - .065 
3 293.1500 5.96416 5.96409 001 5.904244 5.904564 - .005 
3 298.1500 6.38814 6.38329 002 5.904244 5.903725 .009 
3 283.6500 5.16045 5.15947 .019 62 25395K: 6.279469 - .408 
3 285.€500 5.34820 5.34725 .018 6.253355 6.267369 25) 5) 
3 289.1500 5.67198 53712 O11 6.253553 6.25864 -.079 
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. 1500 
- 1500 
. 6500 
.6500 
. 1500 
. 1500 
- 1500 
.6500 
.E50C 
.1500 
. 1500 
. 1500 
.€500 
.6500 
- 1500 
- 1500 
. 1500 
.6500 
- 6500 
- 1500 
. 1500 
- 1500 


283 
285 


K 


1500 
1500 
6500 
6500 
1500 
1500 
1500 


1500 
1500 
1500 
6500 
6500 
1500 
1500 
1500 
6500 
6500 
1500 
1500 
1500 


OIMDVUUNTIVD ADH HIV OHO VU UVIOWUM WYO WU OO WOU ®D OM U1 0101 OD D OF WI UTD D O11 01D D UI G1 U1 O VI i UI O_O GI LF UTM O 


Table II. (Continued) 


PRESSURE 


EXP 


.03833 
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35447 
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.05796 
.52247 
.17651 
.38111 
.73825 
.14315 
.64913 
.17738 
38421 
. 74278 
15210 
.66181 
.18517 
.40368 
.78719 
.22585 
.77549 
. 18830 
.41282 
.80755 
.26115 
.82979 
.19154 
.42238 
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30224 
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.19785 
.44163 
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NOTA 
.00764 
.21258 
.47945 
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.51140 
21383 
.21950 
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.97733 
54496 
.26451 
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.54389 
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30201 
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.21747 
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.464256 
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.6934i 
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.89055 
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Table 


DENS 


. 50000 
-54900 
. 58000 
.60300 
.61300 
.61600 
. 00000 
.95700 
.09900 
. 12500 
. 13600 
. 50000 
.57000 
.61900 
.64800 
.660C00 
. 00000 
.08200 
. 13900 
.171C0 
. 50000 
.59400 
.65900 
.69500 
. C0000 
. 10400 
.17900 
.21900 
. 23300 
. 50000 
.61100 
.69800 
. 74200 
.75700 
.00000 
.11100 
. 21400 
. 26300 
.28100 
. 50000 
.60300 
. 72500 
. 78300 
.80200 
.80900 
-81400 
. 00000 
.08800 
.22700 
. 29800 
.32100 
10.32800 
10.33400 


10.34400 


a od 4 ot 4 
DCODODOHOHHOCHOUOCHKHVOUWUUODDADAAMDAADAANNVINSNINVYVINVANDAADAOAVNANDOANUNANUNUUMS 


298. 
293. 
288. 
284. 
282. 
281 
298. 
293. 
288. 
284. 
282. 
298. 
293. 
288. 
284. 
282. 
298. 
293. 
288. 
284. 
298. 
293. 
288. 
284. 
298. 
293. 
288. 
284. 
282. 
298. 
293. 
288. 
284. 
282 
298. 
293. 
288. 
284. 
282. 
298. 
293. 
288. 
284. 
282. 
281 
281 
298. 
293. 
288. 
284. 
282. 
281. 
281. 


280. 


.650 231.04700 232.02923 -.98223 - 


IV. Comparison of the experimental 
of sound with the formulation. 


JK W. EXP W.CAL DELTA W 


150 241. 91000) | 242596637 -1.05637) == 
150 23259320 OM E2334 20 9 ON mis 2)74eO meas 
150 222.62800 222.53074 .09726 

150 22 N66OO) V2N256N26> 7. 45260 i) 
350 205.51200 207.56345 -2.05145 -. 


.650 201.86200 206.26015 -4.39815 -2. 


150 241.57300 242.43162 -.85862 
150 CSc OOo. Al s(SIl7 2cAtO0Y vot 
150 219.51300 219.41685 .09615 5 
150 205.84400 206.62650 -.78250 -. 
350 194.09100 198.22676 -4.13576 -2. 
150 242.77200 243.43902 -.66702 -. 
150 231.74000 231.73573 -00427 
150 217.77400 217.53422 . 23978 


150 200.18300 200.81873 .-.63573 -. 
350 174.9€900 183.61327 -8.62427 -4. 


150 245.75100 246.27168 -.52068 -. 
150 233.76300 233.61840 - 14460 
150 217.89400 217.48606 .40794 
150 19629800" AI6NS9387, <2 29587) a7 
150 250.88600 251.28200 -.39600 -. 
150 238.05900 237.86640 .19260 
150 220.57500 220.05862 .51638 
150 196.03000 195.83742 .19258 : 
150 258.50200 258.86428 -.36228 -. 
150 245.18700 244.99916 . 18784 
150 226.69700 226.11204 .58496 
150 200.89600 200.24903 - 64697 
350 172.35400 159.84390 12.51010 7 
150 269.04300 269.40842 -.36542 - 
150 255.62100 255.48741 .13359 
150 236.91800 236.41878 -49922 
150 211.68200 210.94564 . 73636 
. 350 185.53700 185.19053 . 34647 
150 282.85400 283.23791 -.38391 - 
150 269.78800 269.55179 . 23621 
150 251.93500 251.43107 -50393 
150 228.79300 228.09131 . 70169 
350 210.08600 210.31108 -.22508 - 
150 300.68200 300.54494 .13706 
150 287.87500 287.18976 .68524 
150 271.83600 271.10114 . 73486 
150 251.84500 250.98933 .85567 
350 238.19300 237.89592 . 29708 


.150 224.90100 229.76267 -4.86167 -2. 
150 322.34300 321.34499 .99801 
150 309.81800 308.21270 1.60530 
150 296.38200 294.77895 1.60305 
150 279.92400 278.32151 1.60249 
350 269.66100 268.41720 1.24380 


650 264.93800 264.00546 93254 
159 261.15200 260.84041 .31159 
150 252.47700 255.96340 -3.48640 -1. 
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speed 


.0766 
- 2580 
- 3220 
. 2584 
.1358 
10523 
SAN Oe 
.3479 
. 1867 
a VGC 
.0876 
. 2000 
. 3067 
Ova 
.0456 
. 2380 
TAOS 
. 3398 
.1247 
.4251 


1617 


. 3096 
.5161 
.5409 
.5725 
-4612 
. 3520 
1193 


3809 
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Cc 
Cc 


aa0n 


MAIN PROGRAM FOR GENERATION OF TABLES IN THIS REPORT. COEFFICIENTS 
ARE OBTAINED FROM BLOCKDATA. 
IMPLICIT. REAL*B (A-H.O-Z) 
REAL TT, ANAMES, QNAMES 
COMMON /COEFS/ A(20).Q(20),ANAMES(20) .QNAMES(20) 
COMMON /OTHER/ COMP ,DPDRHO,DPDTCK 
DIMENSION 1T(27)/276..277.,278.,279. ,280. ,280.5,281.,281.5,282. 
1,282.5,283. ,283.5,284. ,284.5,285. ,286.,287.,283. ,289.,290. ,291. 
2.292 .).293) 294. 296. 296) S007 
159 FORMAT(5X,2HA(,12,2H)=,F12.6,5XA6, 8X2HQ(,12,2H)=,F12.6,5XA6) 
CALL CONGEN 
WRITE (6.159)(I.A(1).ANAMES(1I),.1,Q(1).QNAMES(I),I = 1,20) 
66 READ(5,800,END=99) DENSIN 
B00 FORMAT() 
WRITE(6,44) DENSIN 
44 FORMAT(F12.2,’ MOLE/LITER ISOCHORE‘/111X1H./’ TEMPERATURE PRESSUR 
1E ISOTHERM ISOCHORE INTERNAL ENTHALPY  ENTROPY’SX‘CV’,7X‘CP! 
2, 4X8HVELOCITY/24X DERIVATIVE DERIVATIVE ENERGY’ .41X/OF SOUND’ / 
3 5X'DEG K’,7X’/MPA’ ,5X’L-MPA/MOL MPA/K’,4X'J/MOL’5X5HU/MOL, 
4 4X'U/MOL-K J/MOL-K!’ ,8X‘/M/SEC’/111X1H. ) 
PO 64 II1=1,27 
TEMK=TT (III) 
CALL THERMO(R2KT.R1,7H1,DENSIN, TEMK, PRESS, DVD,CPS,ENTROP,ENTHAL, 
1 CVOL, ENERGY, SOUND) 
DPDRHO=1./DENSIN/COMP 
IF(DABS(TH1).EQ.1.) GO TO 22 
WRITE(6,55) TEMK,PRESS,DPDRHO, DPDTCK, ENERGY, ENTHAL, ENTROP 
1,CVOL,CPS.SOUND : 
55 FORMAT(F10.1,F12.5,F11.4,F10.4,2F10.1,710.2,2F9.2,F9.1,2F10.5) 
GO TO 64 
22 WRITE(6,55) TEMK,PRESS,DPDRHO, DPDTCK, ENERGY, ENTHAL, ENTROP, CVOL 
64 CONTINUE 
GO TO 66 
99 STOP 
END 


ALTERNATIVE MAIN PROGRAM WHICH ENABLES THE USER TO CHOOSE THE SYSTE’m 
OF UNITS FOR THE TABLE. CALCULATION IS ALONG ISOCHORES, WITH THE 
INPUT BEING DENSITY, INITIAL TEMPERATURE, FINAL TEMPERATURE AND 
TEMPERATURE INCREMENT. 

IMPLICIT REAL*B (A-H,O-Z) 
REAL ANAMES, QNAMES 
DOUBLE PRECISION N7,ND,NP.NH 
COMMON /UNITS/ IT,10.1P,1H,NT,ND.NP,NH,FT,FD,FP, FH 
COMMON /COEFS/ A(20),Q(20),ANAMES(20) ,QNAMES(20) 
COMMON /OTHER/ COMP ,DPDRHO, DPDTCK 
159 FORMAT(5X,.2HA(,12,2H)=,F12.6,5XA6,8X2HQ(,12,2H)=,F12.6,5XA6) 
CALL CONGEN 
WRITE (6,159)(1,A(I),ANAMES(I).1,Q(1I),.QNAMES(I),I = 1,20) 
CALL UNIT 
66 WRITE(6,799) 
799 FORMAT(’ INPUT DENSITY, INITIAL T, FINAL T AND T INCREMENT: ’) 
800 FORMATS ) 
READ(S5,800,END=99) DENS,11,72,TI 
WRITE(6,44) DENS,ND,NT.NP,NP,ND,NP,NT,NH,NH,NH,NT,NH,NT 
DENS IN=DENS*FD 
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44 FORMAT(F12.2,2X,A6,’ ISOCHORE’/ 111X1H./’ TEMPERATURE PRESSUR 
1E ISOTHERM ISOCHORE INTERNAL ENTHALPY ENTROPY'5X‘CV‘,7X‘CP 
2‘ .4XBHVELOCITY/24X’DERIVATIVE DERIVATIVE ENERGY’ ,42X’OF SOUND’/ 
3 6X,A6,4X,A6,3X.A4,'/' ,A6,1X,A4,'/',A1,4X,A6,4X,A6, 
BTAX AD we VA tA = 1 All SX, M/SEC L/ 1 411X, 11Hi.) 
TH=T1-TI 
64 TT=TT+TI 
IF(TT.GT.T2) GO TO 66 
TEMK=TTT(TT) ; 
CALL THERMO(R2KT,.R1,TH1,DENSIN, TEMK,PRESS,DVO.CPS,ENTROP,ENTHAL, 
1 CVOL,ENERGY,SOUND) 
DPDRHO=FP/DENS/COMP 
POUT=PRESS*FP 
E=ENERGY* FH 
H=ENTHAL* FH 
S=ENTROP*FH*FT 
CV=CVOL*FH*FT 
CP=CPS* FH*FT 
- DPDT=DPDTCK*FP*FT 
IF(DABS(TH1).EQ.1.) GO TO 22 
WRITE(6,55) TT,POUT,OPDRHO,DPDT,E,H,S,CV.CP,SOUND 
55 FORMAT(F10.2,F12.5,F11.4,F10.4,2F10.1,F11.3,2F9.2,F9.1,2F10.5) 
GO TO 64 
22 WRITE(6,55) TT,POUT,OPDRHO,OPDT,E.H,S,CV 
GO TO 64 
99 STOP 
END 


BLOCKDATA 

DOUBLE PRECISION A,Q 

COMMON /COEFS/A(20) .Q(20) ,ANAMES(20) ,QNAMES(20) 

DATA A/-.007833D0,2*0.0D0.-16.9686D0,5.3439D0, .32500,1.37617D0,0.D0 
1,1.67472D0,23.9809D0.4.82D0, .7662D0, -36.4812D0, -27.35945D0 
2,-12.5003D0,-10.5742D0,4*0.D0/ 

DATA Q/-.008369D0, - .001218D0,28.0541D0, 282. 3452D0,7.634D0,5.0403D0 
1,.1130900,-.777D0.12*0.D0/ 


DATA ANAMES/6HC .6HCOT3 ,6H ,6HCOT2 ,6HDPDT ,6HBETA 
1,6HKO . 6H ,6HB*B .GHA »6HDELTA ,6HK1 » 6EMUC 
2,6HMU1 . 6HMU2 . GHMU3 »6HSOO ,»6HS20 .6HSO1 »6HS21 / 

DATA QNAMES/2*6H »-6HMOL WT,6HTC »6HRHOC  ,6HPC 
1,€HDPCDTC,.6HSLOPDI ,6HP11 »-6HALPHA ,6HPOO «6HP20 ,» 6HP40 
2,6HDELTAI,6HALPHAI ,6HBETA I,6HGAMMAI ,6HPO1 »6HP21 ,6HP41 / 

END 


SUBROUTINE CONGEN 
C THIS SUBROUTINE CALCULATES ALL QUANTITIES NOT DEPENDENT ON R OR THETA 
IMPLICIT DOUBLE PRECISION (A-H,0-Z) 
COMMON /COEFS/ A(20).Q(20) 
DATA DELPP,DELPM.DELPPH, DELPMH/1.D-6,-1.D-6,1.000001D0, .999999D0/ 
EQUIVALENCE (DELI.Q(14)),(DRHOC,A(3)),(BSETA,A(6)),(DELTA,A(11)) 
1,(COTEE2,A(4)),.(DPDT,A(5)),(DELTC,A(8)).(BESQ,A(9)) 
2,(CC,A(1)),(AA,A(10)).(TCK,Q(4)).(RHOC,Q(5)).(PC,Q(6)). (XKO,A(7)) 
3, (ALHI,Q(15)),(BETI.Q(16)),(GAMI,Q(17)).(P01,Q(18)).(P21,.Q(19)) 
4, (ALPHA,Q(10)),(P00,Q(11)),(P20,Q(12)).(P40,Q(13)),(F41,Q(20) ) 
5,(S00,A(17)),(S20,A(18)),(S01,A(19)),(521,A(20) ) 
ALPHA = 2.D0 - A(6)*(A(11) + 1.D0 ) 
IF(DABS(ALPHA).LT.1.D-6) ALPHA=1.D-6 
GAMMA = BETA*(DELTA - 1.D0) 


oy) 


DELI = 0.5000 
IF (DELI.LT.-BETA)DELI=-BETA 

ALHI = ALPHA - DELI 

BETI = BETA + DELI 

GAMI = GAMMA - DELI 
IF(ALHI.LE.DELPPH.AND.ALHI.GE.DELPMH) ALHI=1.000001D0 
IF (ALHI.LE.DELPP.AND.ALHI.GE.DELPM) ALHI = .000001D0 
DELPC = Q(7)*DELTC 
ERR = 2.D0*BETA*DELTA - 1.00 
POO = (BETA*(DELTA-3.D0) -BESQ*ALPHA*GAMMA ) /(2.D0* 

1 BESQ*BESQ*(2.D0-ALPHA)*(1.D0-ALPHA) *ALPHA) 
P20=-(BETA*(DELTA-3. )-BESQ*ALPHA*ERR)/(2.*BESQ*(1.-ALPHA)* (ALPHA) ) 
P40=(ERR-2.)/2./ALPHA 
SOO = (2.-ALPHA)*POO 
S20 = -BETA*(DELTA-3.)/2./BESQ/ALPHA 

IF (XKO.LE.0.D0)XKO = .0001D0 
DA = Q(1) 
DB = Q(2) 
RA = DA/(1.D0 - BESQ) 
RB = DB/(1.D0 - BESQ) 
SWO = 500+S20 
SS1 = AA*XKO*SWO*RA**(1. - ALPHA) 
SS2 = AA*XKO*SWO*RB**(1. - ALPHA) 

D1 = CC*(SS1-SS2) 

PO1 = (BETA*(DELTA-3.D0)-3.DO0*DELI-A(9)*ALHI*GAMI) 

1 /(2.D0*A(9)*A(9)*(2.D0-ALHI)*(1.D0-ALHI)*ALHI ) 

P21 = -((BETA*(DELTA-3.D0)-3.DO*DELI-A(9)*ALHI*ERR) / 

4 (2.D0*A(9)*(1.DO-ALHI) *ALHI) ) 

P41 2 (.5*ERR - 1.)/ALHI 
SO1 = (2.-ALHI)*PO1 
S21 © -(BETA*DELTA-3.*BETI)/2./BESQ/ALHI 
PTW = PO1 + P21 + P41 
SW1 = $01+S21 
SSW1 = AA*A(12)*(RA**(1.DO-ALPHA+DELI) )*SW1 
SSW2 = AA*A(12)*(RB**(1.DO-ALPHA+DELI) )*SW1 
D2 = CC*(SSW1-SSwW2) 
Q(9) = Q(8)*.9905D0 + (D1+D2)/(DB-DA) 
RETURN 
END 


SUBROUTINE THERMO (RHO2KT,R1,TH1,D,7T,PRESS,DVDPT,CP, 

2 ENTROP,ENTHAL,CV, ENERGY , SOUND) 
C CALCULATES ALL THERMODYNAMIC PROPERTIES, BOTH IN 1- AND 2-PHASE 
Cc REGION, FOT GIVEN DENSITY AND TEMPERATURE. 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

COMMON /OTHER/ COMP,DPDRHO,DPDTCK 

COMMON /COEFS/ A(Z20),Q(20) 

DIMENSION S(2),SD(2) 

EQUIVALENCE (COT3,A(2)),(DRHOC,A(3)),(BETA,A(6)),(DELTA,A(11)) 
1,(COT2,A(4)), (DPDT.A(5)), (XMWT,Q(3)), (DELTC,A(8)), (BESQ,A(9Q) ) 
2,(CC,A(1)), (AA, A(10)),(TCK,Q(4)), (RHOC,Q(5)),(PC.Q(6) ) 

3, (ALPHA,Q(10)),(XK1,A(12)), (DELI.Q(14)).(FAC.Q(7)).(XKO,A(7) ) 
4, (ALHI,Q(15)),(P11,.Q(9)),(P1W,Q(18)),(P2W,Q(19)), (P4W,Q(20) ) 
5,(P00,Q(11)), (P20,Q(12)),(P40,Q(13)),(S00,A(17)),(S20,A4(18)) 

IFLAG & O 

XK(1)2A(7) 

XK(2)2A(12) 

8 TER = 2.D0*BETA*DELTA - 1.00 
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TEE=(T-TCK)/TCK 
TSTAR=TEE+1.D0 
DISTIN=1.D0-(1.D0/TSTAR) 
RHO=D/RHOC 
CALL CONVER(RHO,TEE.AMU,TH1.R1,RHO1,S,RHODI,ER) 
TT1=TH1*TH1 
PTW=P1W+P2W*TT14+P4W*TT1*TT1 
P1=AA*XKO*(R1**(2.DO-ALPHA) )*(POO+TT1*(P20+TT1*P40) ) 
P3=RHODI*AMU+PTW*AA*XK1*(R1**(2.DO-ALPHA+DELI) ) 
P4=DTSTIN*(DPDT+DTSTIN*(COT2+DTSTIN*(COT3))) 
PBAR=TSTAR*(P1+P3+P4+1.DO) 
PRESS=PBAR*PC 

IF(IFLAG.EQ.1) RETURN 
PT=DPDT+2.D0*COT2*DTSTIN+3.DO*COT3*DTSTIN*DTSTIN 
UD=A(14)+2.D0*A(15)*DTSTIN+3.DO*A(16)*DTSTIN*DTSTIN 
UFUNCT=S(1)+P11*AMU+PT -RHO*UD+S (2) 
ENERGY=1.D3* (PC*UFUNCT ) /RHO/RHOC 
CALL AUX(R1,TH1,D2PDT2,D2PDMT, RHO2KT, AA,XK,SD,CVCOEX) 
IF(DABS(TH1).EQ.1.D0) RHO2KT=1.D30 
DVDPT=- (( RHO2KT *RHOC**2) /((D**3)*PC*TSTAR) ) 
COMP =RHO2KT/PC/RHO* *2 
DPDRHO=1./D/COMP 
CV1=2.DO*P11*D2PDMT/RHO2KT 
CV2= (D2PDMT*D2PDMT ) /RHO2KT 
D2PZER=2.D0*COT2+6*COT3*DTSTIN 
D2MZER=2.D0*A(15)+6.DC0*A(16)*DTSTIN 
CV3=RHO*D2MZER 
CV4=P11*P11/RHC2KT 
CVS= (D2PZER-CV3-CV1-CV2-CV4+D2PDT2)/TSTAR**2 
CVCOEX=(CVCOEX*AA) /((1.00-BESQ) **2) 
IF(DABS(TH1).EQ.1.D0) CVS=(D2PZER-CV3+CVCOEX) /TSTAR**2 
CV=1.D3*PC*CVS/D/TCK 
TAU=(-1.D0)/TSTAR 
DPDTR=PT+P11* ( AMU- ( RHO/RHO2KT) )+S(1)+S(2) -RHO*D2PDMT/RHO2KT 
DPDTR=PBAR- (TAU*DPDTR) 
DPDTCK=DPDTR*PC/TCK 
CP=(RHO2KT*DPDTR**2/RHO**2+CVS) *PC*1.D3/D/TCK 
ENHOL=PRESS/PC/TSTAR- ( TAU*UFUNCT ) 
ENTHAL=1.D3*PC*TSTAR* ENHOL/RHOC/RHO 
UZERO=A(13)+DTSTIN*(A(14)+DTSTIN*(A(15)+DTSTIN*A(16))) 
IF(DABS(TH1).EQ.1.C0) GO TO 44 
SOUND=-CP/CV/D/D/DVDPT/XMWT 
SOUND=DSQRT (DABS(SOUND) )*1.D3 
ENTROP=(ENHOL-RHO* (UZERO+AMU) ) *PC*1.D3/RHO/RHOC/TCK 
CONT INUE 
RETURN 
ENTRY PRES(R1,TH1,D.7T. PRESS) 
IFLAG = 1 
GO TO 8 
END 
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SUBROUTINE AUX(R1,TH1,D02PDT2,D2PDMT ,D2PDM2,AA,xXK,SD,CVCOEX) 


C CALCULATES SECOND DERIVATIVES. IS CALLED BY THERMO 


30 


IMPLICIT REAL*8 (A-H,O-Z) 
COMMON /COEFS/ A(20),Q(20) 
DIMENSION XK(2),.S(2), SD(2),W(2),.Y¥(2),2(2) ,COEX(2) 
EQUIVALENCE (CC.A(1)),(BETA,.A(6)).(BESQ,A(9)), (DELTA, A(11)) 
1, (ALPHA,Q(10)),(S00,A(17)),(S20,A(18)).(S01,A(19)),(S21,A(20) ) 
DELI = 0O.DO 
§(1)=S00+S20*TH1*TH1 
$(2)=S01+S21*TH1*TH1 


SD(1) = 2.*TH1*S20 
SD(2) = 2.*TH1*S21 
Ww = 0.DO 
YY = 0.DO 
ZZ = Q.DO 
GAMMA = BETA*(CELTA - 1.D0) 
TT1 = TH1*TH1 


TER = 2.DO*BETA*DELTA - 1.D0 
G = (1.+(BESQ*TER-3.)*TT1 - BESQ*(TER-2.)*TT1*TT1) 
CVCOEX = 0.DO 

DOMNSO; I=) 2 


ALHI = ALPHA - DELI 

BETI = BETA + DELI 

GAMI = GAMMA - DELI 
W(1)=(1.-ALHI)*(1.-3.*TT1)*S(1)-BETA*DELTA*(1.-TT1)*THI*SD(I) 
W(1) (W(1) *(R1**(-ALHI)))/G 


W(I) = W(1)*XK(1) 
WW = WW + W(I) { 
Y(I) = BETI*(1.D00-3.D0*TT1)*TH1 - BETA*DELTA*(1.D00-TT1)*TH1 
) (Y(I) *(R1**(BETI - 1.D0)))*XK(1)/G 
= Vee SA(ES)) 
Z(I) = 1.D0 - BESQ*(1.D0 - (2.DO*BETI))*TT1 
) 


<u 


= (Z(I) *(R1**(-GAMI)))*XK(1)/G 
ZL) HZ(Oln) 
Al (BETA*(DELTA-3.D0)-3.D0*DELI-BESQ*ALHI*GAMI ) 
1 /(2.00*BESQ*BESQ*(2.D0-ALHI)*(1.DO-ALHI)*ALHI ) 
A2 = 1+((BETA*(DELTA-3.D0)-3.DO*DELI-BESQ*ALHI*TER)/ 
1 (2.D0*BESQ*(1.D0-ALHI)*ALHI) ) 
3 cor) / \V2) 
A4 = 1.D0+((TER-2.D0)/(2.D0*ALHI ) ) 
Fi = Ai + A2:+ AG 
COEX(I) = ((2.00 - ALHI)*(1.D0 - ALHI)*(R1**(-ALHI) )*FI*XK(T)) 
CVCOEX = CVCOEX + COEX(1) 
DELI = 0.5D0 


CONTINUE 
D2PDT2 = AA*WW 
D2POMT = YY + AA*CC*WW 
D2PDM2 = ZZ/AA + 2.DO*CC*YY + (CC#*2)*AA*WW 
RETURN 
END 


SUBROUTINE SS(R,TH,S,SD) 


C CALCULATES THE PARAMETRIC PART OF THE ENTROPY FUNCTION 


Cc 


FOR GIVEN R AND THETA. 


IMPLICIT DOUBLE PRECISION(A-H,O-Z) 

DIMENSION $(2).SD(2) 

COMMON /COEFS/ A(20),Q(20) 

EQUIVALENCE (ALPHA,Q(10)),(BETA,A(6)),(B8ESQ,A(9)),(DELTA,A(11)) 


1,(0ELI,Q(14)), (ALHI.,Q(15)),(BETI,Q(16)).(GAMI,Q(17)),(P00,Q(11)) 
2,(P01,Q(18)),(S00,A(17)),(S20.A(18)),(S01,A(19)),(S21,A(20) ) 


TT = TH*TH 

S$(1)= SOO + S20*TT 

SD(1) = 2.*S20*TH 

S(2) = SO1 + $21*TT 

SD(2) = 2.*S21*TH 
S(1)=S(1)*A(10)*A(7)*R**(1.-ALPHA) 
$(2)=S(2)*A(10)*A(12)*R**(1.-ALHI) 
RETURN 

END 
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SUBROUTINE CONVER(RHO,TEE,AMU,TH1,R1,RHO1S,S1,RHODI,ERROR1) 
IMPLICIT REAL*8 (A-H,0-Z) 
COMMON /COEFS/ A(20),Q(20) 
DIMENSION $1(2),SD(2) 
EQUIVALENCE(BETA,A(6)),(DELTA,A(11)),(XK1,A(12)),(CC,A(1)) 
1 ,(ALHI,Q(15)), (ALPHA,Q(10)),(BESQ,A(9)),(P11,.Q(9)) 
2 ,(DELI,Q(14)),(P1W,C(18) ).(P2W,Q(19)), (P4W,Q( 20) ) 
3,(AA.A(10)),(XKO,A(7)),(S00,A(17)),(S20,A(18)),(BETI,Q(16) ) 
TSTAR = TEE+1.D0 


DTSTIN = 1.D0 - (1.DO/TSTAR) 
R1=DTSTIN 5 
IF(DTSTIN.LT.O.) R1=DTSTIN/(1.-BESQ) 
TH1=0. 


MR CDMSTEUN SST Ore)! nite. 
CALE PSSICRAT Ty S|) SiD)) 
RHODI = 1.D0 + P11*D7STIN 
RHODIT = RHODI + CC*(S1(1)+S1(2)) 
DRHO = RHO - RHODIT 
AMU = 0.DO 
IF(DTSTIN.GT.O.DO) GO TO 1 
RHO1CO = XKO*R1**BETA + XK1*R1**BETI 
TWOFAZ = RHO1CO 
IF (DABS(DRHO).GT.TWOFAZ) GO TO 1 
RHO1S = DSIGN(RHO1CO,DRHO) + CC*S1(1) 
TH1 = DSIGN(1.D0,DRHO) 
ERROR1 = 1.D0 
GO TO 999 
1 CONTINUE 
IF(DRHO.NE.O.DO) GO TO 2 
TH1 = O.DO 
R1 = DTSTIN 
RHO1S = CC*S1(1) 
2 CONTINUE 
C RULE FOR FIRST PASS 
Y1 = DTSTIN 

DEN1 = RHO - RHODIT 
CALL RTHETA(R1,.1H1,DEN1,Y1) 
TT=TH1*TH1 

AMU = AA* R1**(BETA*DELTA)*TH1*(1.90-TT) 
Y1 = DTSTIN + CC*AMU 
CABERSSIGRNS Tht Silt Si) 

RHOWEG = XK1*(R1**(DELI+BETA))*TH1 + CC*S1(2) 
RHO1S = DEN1 + CC*S1(1) + RHOWEG 
ERROR1 = RHO - RHODI - RHO15 
IF( DABS(ERROR1).LT.1.D-5 ) GO TO 999 

C RULE FOR SECOND PASS 


DEN12 = RHO - RHODI - CC*S1(1) + RHOWEG 

IF (DEN12.EQ.DEN1) DEN12 = DENI - 1.0-06 
CALL RTHETA(R1,TH1,DEN12,Y1) 
TT = TH1*TH1 

AMU = AA* R1**(BETA*DELTA) *TH1*(1.D0-TT) 
Y1 = DTSTIN + CC*AMU 
CALL SS(R1,TH1,S1,SD) 

RHOWEG = XK1*R1**(DELI+BETA)*TH1 + CC*S1(2) 
RHO1S2 = DEN12 + CC*S1(1) + RHOWEG 
ERROR2 = RHO - RHODI - RHO1S2 
IF( DABS(ERROR2).GT.1.0-5 ) GO TO 998 
ERROR1 = ERROR2 
RHO1S = RHO1S2 
GO TO 999 

998 CONTINUE 
C RULE FOR NTH PASS 

DEN2 = DENi2 
DO 44 ISIG = 1,10 
SLOPE = (ERROR2-ERROR1)/(DEN2 -DEN1) 

HOLD = DEN2 
DEN2 = DEN1 - (ERROR1/SLOPE) 
CALL RTHETA(R1,TH1,DEN2.Y1) 
TT = THI*THI 
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44 


66 
999 


cc 
cc 
cc 


100 


500 
601 
498 


495 


496 


609 


450 


AMU = AA* R1**(BETA*DELTA)*TH1*(1.D0-TT) 

Y1 = DTSTIN + CC*AMU 

CALL SS(R1,TH1,S1,SD) 

RHOWEG = XK1*R1**(DELI+BETA)*TH1 + CC*S1(2) 
RHO1S = DEN2 + CC*S1(1) + RHOWEG 

ERROR1 = ERROR2 

ERROR2 = RHO - RHODI - RHO1S 

IF( DABS(ERROR2).LT.1.D-6 ) GO TO 999 

DEN1 = HOLD 
CONTINUE 

IF(DABS(ERROR2).GT.1.D-6) WRITE (6,66) ERROR2,RHO,DTSTIN,DEN2,CC 
FORMAT(1X,8F9.5) 
CONTINUE 

RETURN 

END 


SUBROUTINE RTHETA(R,THETA, RHO,TEE ) 

THIS VERSION WILL FIT DATA FOR THETA GT 1.0 AND LT 1.000001 AS IF M 
SOLVES RHO = EM*THETA*(R**BETA) 
AND TEE = R*(1.D0-BESQ*THETA*THETA) 

IMPLICIT REAL*8 (A-H,O-Z) 

COMMON /COEFS/ A(20).Q(20),QQ(20) 

EQUIVALENCE (BETA,A(6)),(EM,A(7)),.(BESQ,.A(9)) 
IF(EM.LE.O.DO .OR. BESQ.LE.1.00) GO TO 600 

ABSRHO = DABS( RHO ) 

IF(A8SRHO .LT. 1.D-12 ) GQ TO 690 

BEE = DSQRT(BESQ) 

IF(DABS(TEE) .LT. 1.D-12) GO TO 495 

IF( TEE .LT. O.DO ) Z = 1.00-(1.D0-BEE)*TEE/(1.D0-BESQ) 
1 *(EM/ABSRHO)**(1.D0/BETA) 

IF(TEE.GT.O.DO0) Z = (1.D0+TEE*(EM/BEE/ABSRHO)**(1.DO0/BETA))**-BETA 
IF( Z.GT.1.00234D0*BEE ) GO TO 496 

C = -RHO*BEE/EM/DABS(TEE)**BETA 

Z = DSIGN(Z,RHO) 

DO 500 N = 1, 16 

Z2 = FA4Z 
ZZ = 1 IDO Z2 


DZ Z3*(Z+C*DABS(Z3)**BETA)/(Z3+2.D0*BETA*Z2) 
PAE Low [DS 

IF(DABS(DZ/Z) .LT. 1.D-12) GC TO 498 

CONTINUE 


WRITE(6,450) R,THETA, RHO,TEE,BETA,EM,BESQ 
IF(DABS(THETA).GT.1.0001) THETA=THETA/DABS( THETA) 
RETURN 
THETA = Z/BEE 
R = TEE/(1.D0-Z*Z) 
R = DABS(R) 
RETURN 


THETA = DSIGN(1.DO,RHO)/BEE 
R = (RHO/(EM*THETA) )**(1.D0/BETA) 
RETURN 
THETA = DSIGN(1..RHO) 
R = TEE / ( 1.000 - BESQ ) 
R = DABS(R) 
RETURN 


IF(DABS(TEE) .LT. 1.D-12) GO TO 60) 
IF(TEE.LT.O.DO) GO TO 496 
THETA = 1.D-12 
R = TEE 
FORMAT(2X,’RTHETA DOES NOT CONVERGE’ , 7(1X, E12.6) ) 
RETURN 
END 
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Cc 
Cc 


SUBROUTINE SATDEN CALCULATES THE VAPOR PRESSURE AND THE 
SATURATED LIQUID AND VAPOR DENSITIES AT THE INPUT TEMPERATURE 
SUBROUTINE SATDEN(T,DL,DV,PRESS) 
IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

COMMON /OTHER/ COMP,CPDRHO,DPDTCK 

COMMON /COEFS/ A(20).Q(20) 

DIMENSION S(2).SD(2) 

EQUIVALENCE (COT3,4(2)),(ORHOC,A(3)),(BETA,A(6)),(DELTA,A(11)) 
1,(COT2,A(4)),(DPDT,A(S)), (XMWT,Q(3)), (DELTC,A(8)),(BESQ,A(9) ) 
2,(CC,A(1)).(AA,.A(10)).(TCK,Q(4)),(RHOC,Q(5)),(PC,Q(6)) 

3, (ALPHA,Q(10)),(XK1,A(12)), (DELI,.Q(14)).(FAC,Q(7)). (XKO,A(7)) 
4,(ALHI,Q(15)),(P11.0(9)).,(P1W,Q(18)). (P2W.Q(19)),(P4W,Q(20)) 
5,(P00,Q(11)),(P20,.Q0(12)),(P40,Q(13)),(S00.A(17)),(S20,A(18) ) 
6,(S01,A4(19)),(S21,.A(20)), (BETI,Q(16)) 

ERR = 2.DO*BETA*DELTA - 1.D0 

TEE=(T-TCK)/TCK 

TSTAR=TEE+1.DO 

DISTIN=1.D0-(1.DO/TSTAR) 

Di=1.+P11*DTSTIN 

RR=DTSTIN/(1.-BESQ) 

D2=XKO*RR**BETA + XK1*RR**BETI 

D3S=CC*AA* (XKO*(S00+S20)*RR**(1.-ALPHA) + XK1*(SO14+S21) * 

1 RR**(1.-ALHI)) 

RHOL=D1+D2+D3 

RHOV=D1-D2+D3 

DL=RHOL*RHOC 

DV=RHOV*RHOC 

CALL CONVER(RHOL,TEE,.AMU,TH1,R1,RHO1,S,RHODI,ER) 

PTW=P1W+P2W+P4W 

P1=AA*XKO*(R1**(2.00-ALPHA) )*(PO00+P20+P40) 

P3=RHODI*AMU+PTW*AAYXK1*(R1**(2.D0-ALPHA+DELI ) ) 

P4=DTSTIN*(DPDT+DTSTIN*(COT2+DTSTIN*(COT3))) 

PBAR=TSTAR* (P14+23+P4+1.D0) 

PRESS=PBAR*PC 

RETURN 

END 


SUBROUTINE DFIND(D,P,DGUESS,T,DVDPT) 


C THIS SUBROUTINE WILL CALCULATE BY ITERATION THE DENSITY FOR A 
Cc GIVEN PRESSURE AND TEMPERATURE. AN INITIAL GUESS FOR THE 
Cc DENSITY ("DGUESS") IS ALSO REQUIRED. 
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IMPLICIT DOUBLE PRECISI1ON(A-H,0-Z) 
DIMENSION S(2),SD(2),XK(2) 

COMMON /COEFS/ A(20).Q(10) 
XK(1)=A(7) 

XK(2)=A(12) 

DD = DGUESS 

VOM t= DO 

Wu So 

DYLIM = .0001 

DOnOmuate. 2 

CALL PRES(R1,TH1,DD,TT, PRESS) 
DY=(P-PRESS) 

CALL AUX(R1,TH1,D2PDT2,D2PDMT, RHO2KT,A(10),XK,S,SD,CVCOEX) 
TSTAR=T/(Q(4)+A(8)) 

DVDP=- ((RHO2KT*(Q(5)+A(3))**2)/( (DD**3)*(Q(6)+O(7)*A(8))*TSTAR) ) 
IF(DVDP.LT.-2.00) DVDP=-2.D0 
IF(DVDP.GT.1.D-2) DVDP=1.0-2 

DELV = DY*DVDP 

VOL = VOLI + DELV*.6D0 

DD = 1.DO/VOL 

VOLI = 1.DO /DD 

IF(DABS(DY).LT.DYLIM) GO TO 81 
CONTINUE 
WRITE(6,1) P,PRESS,DGUESS,DD,DVOP.DELV 
FORMAT(‘’ NO CONVERGENCE: ’ ,6F12.5) 

D = DD 
RETURN 
END 
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TH 


TH 
I 
A 
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aaan (one) 


30 


31 


32 


33 


99 
10 
11 
12 
13 
14 
15 


C FU 


E FOLLOWING 2 SUBROUTINES FORM A PACKAGE ALLOWING THE USER TO 
OPERATE IN A SYSTEM OF UNITS OF HIS CHOICE. 
SUBROUTINE UNIT 


IS SUBROUTINE QUERIES THE USER AS TO HIS CHOICE OF UNITS AND SETS 
NTERNAL PARAMETERS APPROPRIATELY. THE INTERNAL UNITS OF THE PROGRAM 
RE TEMPERATURES IN K. DENSITIES IN MOL/DM3, ALL OTHER QUANTITIES 
ALCULATED IN DIMENSIONLESS UNITS AND DIMENSIONED AT OUTPUT TIME. 


IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

DOUBLE PRECISION NT,ND,NP,NH,NNT,NNO,NNP,NNH 

COMMON /UNITS/ 1T,1D.1P,1H,NT,ND,NP,NH,FT,FO,FP,FH 
DIMENSION FFD(4),FFP(5),FFH(6),NNT(4) ,NND(4),NNP(5).NNH(6) 
DATA FFD/.03564541D0,35.6454100,1.D0, .5709853D0/ 

DATA FFP/1.00,10.D0,9.869232667D0,145.037738D0,10.1971D0/ 


DATA FFH/2* .03564541D0,1.D0,8.51379D-3, .2388459D0, .01532476D0/ 


DATA NNT/1HK,1HKC,1HR, 1HF/ 
DATA NND/6HKG/M3 .6HG/CM3 ,6HMOL/L ,6HLB/FT3/ 
DATA NNP/6H MPA ,6H BAR’ ,6H ATM’ ,6H PSI ,6HKG/CM2/ 


DATA NNH/6HKU/KG .6H J/G- ,6HJ/MOL ,6HCAL/G ,7HCAL/MOL,6HBTU/LB/ 


DATA A1,A2,A3,A4/8HTEMPERAT ,7HDENSITY 
1,8HPRESSURE, SHENERGY / 
WRITE(6,11) Ai 

WRITE(6,12) 

READ(5,10,END=99) IT 

IF(IT.EQ.0) STOP 

IF‘1T.GT.4) GOTO 30 

NT=NNT(IT) 

WRITE(6,11) A2 

WRITE(6,13) 

READ(5,10,END=99) ID 

IF(ID.GT.4 .OR. ID.LT.1) GOTO 31 
ND=NND( 1D) 

FD=FFO( ID) 

WRITE(6,11) AZ 

WRITE(6,14) 

READ(S5,10,END=99) IP 

IF(IP.GT.5 .OR. IP.LT.1) GOTO 32 
NP=NNP( IP) 

FP=FFP(IP) 

WRITE(6,11). A4 

WRITE(S,15) 

READ(5,10,END=99) IH 

IF(IH.GT.6 .OR. IH.LT.1) GOTO 33 
NH=NNH( IH) 

FH=FFH( IH) 

RETURN 

STOP 

FORMAT ( ) 

. FORMAT(’ ENTER UNITS CHOSEN FOR ’,A8) 
FORMAT(’ CHOOSE FROM 1=DEG K, 2=DEG C, 3=DEG R, 4=DEG F’) 
FORMAT(’ CHOOSE FROM 1=KG/M3, 2=G/CM3, 3=MOL/L. 4=LB/FT3' ) 
FORMAT(’ CHOOSE FROM 1=MPA, 2=BAR, 3=ATM, 4=PSIA, 5=KG/CM2’) 


FORMAT(’ CHOOSE FROM 1=KU/KG, 2=JU/G, 3=J/MOL, 4=CALORIES/G, 5=CALO 


{RIES/MOL, 6=BTU/LB’ ) 
END. 
FUNCTION TTT(T) 
NCTION TO CONVERT INPUT TEMPERATURES IN EXTERNAL UNITS TO DEG K 
DOUBLE PRECISION T,TTT,FT,FO,FP,FH,NT,NO,NP,NH 
COMMON /UNITS/ I17,1D,1P,1H,NT,NO,NP,NH,FT,FO,FP,FH 
GO T0 .(1.2,3,4).1T ; 
VU VS ; 
Falc=ate 
RETURN 
TTT=T+273.15D0 
Fane 
RETURN 
TTT=T/1.8D0 
FT=.5555555555556D0 
RETURN 
TTT=(T+459.6700)/1.8D0 
FT=.5555555555556D0 
RETURN 
END 66 


Appendix E. Units and Conversion Factors 


Units used in the tables 


Pressure - MPa 
Temperature - K 

Density - mol/dm? 

Energy, Enthalpy - J/mol 
Entropy - J/(mol-K) 
Specific heat - J/(mol-K) 
Speed of sound - m/s 


28.0541 


Molecular weight of ethylene M 
R = 8.31434 J/(mol 


Gas constant* 


Reduction Factors 


Pressure 

MPa to 
Pa 106 
bar 10 
atmosphere 9.869233 
psi, Ibf/in2 145.037738 

Density 

mol/dm? to 
kg/m? 28.0541 
Tbm/ft? 1. 7513585 


Energy, Enthalpy 

J/mol to Multiply table value b 
kJ/kg | 0.0356454 
BTU/1bm 0.01532476 


Entropy, Specific heats 
J/(mo1l.K) to Multiply table value b 


kJ/(kg.k) 0.0356454 
BTU/(1bm. R) 0.00851375 
Temperature 
K to Transform table value b 


G GeO) =. (Khe 273.15 
a (U5) 2 (0). ecole 2 
R (T.R) = g (TK) 


Speed of sound 


m/s to Multiply table value by 
ft/s 3.28083 


*The actual value of the gas constant is not required for the programs in this 
Tech Note. The value given here is that used in reducing the original Burnett 
PVT data (ref. 3, 18) and in NBS Tech Note 1045 (ref. 5). 
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Temperature, K 
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Fig. 1 Range of validity of the present formulation. Left: in 
temperature - density space. Right: in pressure-temperature 
space. : 
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Density Deviations 
T = 298.15 K © 
4 DH 
x TWW 
y+ HLSB 
o TZ 
MJ 


Ap, % 


Density, mol/dm* 


PuGe 2 Density deviation, in %, between the experimental PVT data 
and the present formulation, if compared at the same pressure 
at 298.15 K. <A DH [10]; x TWW [11]; + HLSB [3], 99.993% pure; 
Y HLSB [3], 99.999% pure; o TZ [12,13]; —— Md [5]. 


Pressure Deviations 
T = 298.15 K 


a ee 
eta lacs : 
7 x 8¢e 11 


TWW 


Rig. 3 Pressure deviations, in %, between the experimental PVT data 
and the present formulation, if compared at the same density 
at 298.15 K. Symbols as in Fig. 2 


69 


4 Density Deviations 


T = 293.15 K 
2 
Density, mol/dm? 
0 ; : cid ian S 
se YOO um. oy 7 cs 8° J a 10 11 
=~ A r - 
<4 p A 
—.2 
4 DH 
TWW 
4 y+ HLSB 
Oe WZ 
MJ 
Fig. 4 As in Fig. 2, but at 293.15 K. 
ti Pressure Deviations 
T = 293.15 K 
1 
0 pACte ts 
: 6 7 
on 
a 
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—.| 
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_9 y+ HLSB 
o7 TZ 
MJ 


Fig. 5 As in Fig. 3, but at 293.15 K. 
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Density Deviations 


. Ob 287<T< 290 K 
« TWW 
r+ HLSB 
ou TZ 


Density, mol/dm° 


Fig. 6 As in Fig. 2, but for temperatures from 287-290 K. 


Pressure Deviations 
287< T< 290 K 


| 
xe A 
2. 
= + 
0 


Density, mol/dm* 


Fig. 7 As in Fig. 3, but for temperatures from 287-290 K. 
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25 Density Deviations 
284<T<287 K 


Fig. 8 As in Fig. 2, but for temperatures from 284-287 K. 


Pressure Deviations 
284<T< 287 K 


Density, mol/dm? 


Fig. 9 As in Fig. 3, but for temperatures from 284-287 K. 
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Fig. 10 
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Density Deviations 
278<T<284 K 


Density, mol/dm* 


4DH 

x TWW 
v + HLSB 

o7Z 


As in Fig. 2, but for temperatures from 278-284 K, Uc@os 
spanning the critical temperature. The large density departures 
between the different data sets are caused by the divergence 

of the compressibility and the expansion coefficient. 


Pressure Deviations 
218<T<284 K 2 DH 


x TWW 
Y + HLSB 
° TZ 


6.0 


—.02 


Fig. 11 


Density, mol/dm? 


Pressure deviations, in %, between experimental PVT data and 
the present formulation, if compared at the same density. The 
comparison is for the subcritical densities, 5.5-6.5 mol/dm3, 
and for the temperatures from 278-284 K. Symbols as in Fig. 2. 
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Ap, % 


04 


.04 


Rallgeals 


2 DH Pressure Deviations 
x TWW 218 <1T< 284 K 
Y + HLSB 


° 3 
Density, mol/dm 
Fig. 12 As in Fig. 11, but for densities from 6.5-8.5 mol/dm?, i.e., 
spanning the critical density. 


Pressure Deviations 
218<T<284 K 


Density, mol/dm? x : 


As in Fig. 12, but for liquid-like densities from 8.5 - 10.5 
mol/dm?. 


002 VAPOR PRESSURE DEVIATIONS 
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Fig. 14 Departures of the vapor pressures from the present formulation 
x DH [10]; @ HLSB [3], 99.993% pure sample; a HLSB [3], 
99.999% pure sample. 
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Fig. 15 Departures of the coexisting densities from the present formulation. 
Upward-pointing line segments denote vapor, downward-pointing 


ones liquid densities. 
&® DH [10]; 0 HLSB, 99.993% pure sample [3], o HLSB, 99.999% 


pure sample [3]. 
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Fig. 16 The speed of sound as measured by Gammon [4], 0, compared with 
full curves. Only below Ty very near the phase boundary do 


systematic differences exist. 
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Fig. 17. Departures of the speed of sound data [4], in %, from the 
prediction of our surface. Some points have very large departures; 
their estimated error, however, is also quite large. Near the 
phase boundary, however, the departures of 1-2% are systematic 
and significant. 
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The difference between the prediction of McCarty and Jacobsen 
for several thermodynamic properties and the present formulation. 
Top: enthalpy difference AH; center: energy difference AS. 

The M.J. formulation has a somewhat different density dependence 
than the present one. Substantial differences occur near the 
critical density and within a few K from the critical temperature. 


Fig. 


78 


+20% ACs 


283.15 K 


+10% 


Fig. 20 The percentage difference between second-derivative values 
according to the McCarty-Jacobsen formulation and the present 
one. Top: speed of sound difference AC_; center: difference 
in specific heat at constant volume AC; bottom: difference 
in specific heat at constant pressure AC.. The differences 
between the two formulations are substantial within the range 
of 5 K from the critical point. The analytic formulation M.J. 
will have a non-zero thermodynamic speed of sound and a finite 
C, at the critical point, while it underestimates the strength 
of the divergence of C: 
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Fig. 21. The region of P-T space that should be avoided in custody 
transfer if the pressure is measured to, respectively, 0.1%, 
0.05% or 0.02% and an accuracy of 0.1% is desired in density. 
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transfer if the temperature is measured to, respectively, 


0.1%, 0.05% and 0.02% and an accuracy of 0.1% is desired in 
density. 
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Fig. 23 The region of P-T space that should be avoided in custody 
respectively, 0.005%, 0.01% or 0.02% (in mole fraction) of an 


impurity very different in volatility from ethylene, is present, 
and an accuracy of 0.1% is desired in density. 
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